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ABSTRACT 


An overview and investigation of the more popular digital filter design tech- 
niques are presented, with the intent of providing the filter design engineer a com- 
plete and concise source of information. Advantages and disadvantages of the 
various techniques are discussed, and extensive design examples used to illustrate 
their application to specific design problems. Both IIR (Butterworth, Chebyshev 
and elliptic), and FIR (Fourier coefficient design, windows and frequency sampling) 
design methods are featured, as well as the Optimum FIR Filter Design Program 
of Parks and McClellan, and the Minimum p - Error IIR Filter Design Method of 


Deczky. 
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I. INTRODUCTION 


Digital filter design today is a rapidly growing field with myriad applications 
in the areas of signal processing, image processing, filtering, prediction, and estima- 
tion. Extensive research has resulted in a plethora of available information, which 
can be overwhelming to the engineer presented with a specific design problem. As 
the title suggests, the purpose of this thesis is to present an overview of the more 
popular design techniques in an attempt to consolidate the information available in 
the literature, and provide an easily readable reference manual. Detailed examples 
are given to illustrate the application of selected techniques to specific filter design 
problems. 

Due to the wealth of information available, a summary and investigation of 
every filter design technique is not possible. As stated the more popular meth- 
ods are outlined in detail, however, references are included as sources of further 
information. 

The survey is threefo'd, consisting first of recursive IIR design methods, fol- 
lowed by nonrecursive FIR design, and finally computer-aided design (CAD). A 
summary of the chapter contents follows to enable the reader to immediately lo- 
cate the section applicable to his/her particular design problem. 

Chapter II, Recursive Filter Design, presents Butterworth, Chebyshev and 
elliptic filter design from both a traditional approach, wherein analog prototype 
filters based on design specifications are converted to digital versions using the 
bilinear transformation, and a direct design approach that eliminates the need to 


determine an analog prototype. 


12 


Chapter III, Nonrecursive Filter Design, includes an overview of Fourier co- 
efficient design, windowing, and a derivation and illustration of the method of 
frequency sampling. Among the design examples presented is the design of a band- 
pass differentiator and a bandpass integrator. 

In Chapter IV, Computer-Aided Design, the Optimum FIR Filter Design 
Method of Parks and McClellan (that employs the Remez Exchange Algorithm) 
is investigated and applied to a bandpass filter design problem. Also included is 
a short discussion of a method to enhance the application of this program to the 
design of high-order filters. 

To complete the chapter, the Minimum p-Error Design Method for IIR Filter 
Design, that utilizes the Fletcher-Powell algorithm, is presented. Furthermore, a 
discussion of the advantages and disadvantages of these two iterative techniques is 


included. 
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II. RECURSIVE FILTER DESIGN 


A. INTRODUCTION 

The recursive realization of digital filters is advantageous, in that the desired 
frequency response can be obtained using a lower order filter than if a nonrecursive 
realization were used (assuming linear phase is not required). This is because the 
filter frequency response is influenced by both the poles and zeros of the filter 
transfer function, whereas the frequency response of a nonrecursive realization is 
determined only by the filter’s zeros. 

Traditionally, the design of recursive digital filters involves the determination 
of an analog prototype using one of several methods: Butterworth, Chebyshev or 
elliptic, to name a few, and converting these to a digital version using impulse- 
invariant design or bilinear transformation methods, the latter being the more 
popular. 

This chapter first presents a review of the traditional design techniques us- 
ing examples involving Butterworth, Chebyshev and elliptic filters; then a Direct 
Design method is presented, whereby the requirement for determining an analog 
prototype is eliminated, thus reducing the algebraic complexity of recursive filter 


design. 


B. TRADITIONAL DESIGN TECHNIQUES 

Recursive filter design techniques, as stated in the introduction, involve deter- 
mining an analog prototype filter from given filter specifications, and then convert- 
ing the prototype to a digital version using a bilinear transformation. Since this is 


a well known design procedure [1], details of the derivation will not be considered. 
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However, as a quick review, three examples involving the design of Butterworth, 
Chebyshev and elliptic filters will be given. 

It should be noted that the tables of analog prototype transfer functions found 
in this section (Tables 2.2, 2.3 and 2.4), all have a critical frequency of unity 
(w, =1). Thus, to obtain a filter transfer function with a different critical fre- 
quency based on these prototypes requires the use of an appropriate analog fre- 
quency transformation (Table 2.1). 

For example, suppose a lowpass analog filter with a critical frequency of w, 


is desired. The following relationship applies: 


Hyp (jw,) = Hxpp(j1) (2.1) 


where Hyp, is the prototype filter with a critical frequency of w, = 1. 
To convert the prototype transfer function to the desired transfer function, 
Hy, p(s), the s in the lowpass prototype filter is replaced by s/w, as: 
Hy p(s) = Hzpp(s) (2.2) 


= 
I= 
Cc 


Table 2.1 (Reference 1) lists the appropriate frequency transformations for high- 
pass, bandpass and bandstop filters. A summary of the design procedure follows. 
1. Bilinear Transformation Design 

From given analog filter specifications, determine the approppriate s- 
domain design technique (i.e., Butterworth, Chebyshev or elliptic). 

Design the analog prototype filter in accordance with the technique se- 
lected above. This involves the translation of the filter specifications to those of 
a lowpass prototype. The lowpass prototype is then designed according to the 


translated specifications. 
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Transform the lowpass prototype to the desired lowpass, highpass, band- 
pass, or bandstop filter according to Table 2.1. 
Transform the analog filter transfer function to a digital version using the 


bilinear transformation. 


H(z)\= H(s) (233) 
g= i=l 
z+1 
Example 2.1 Butterworth Filter 
Design a digital filter that is flat in the passband from 0 to the 3 dB cutoff 
frequency, f,, of 2 kHz. For frequencies greater than 4 kHz, the attenuation should 
be at least 10 dB. The sampling rxtc is 20 kHz. 
Step 1: A Butterworth design is called for because a flat passband is desired. 


Step 2: Convert the critical analog design frequencies to digital. 
Sampling frequency: f,=20kHz => w, = 40x 10°m rad/s 
Cutoff frequency: fo =2 kHz = Sw — 4 10 aad 


Stopband frequency: f, =4kHz => w, =8 x 10° rad/s 


(4 x 10°) 
O. =) — ee = “90 x 103 = eo rad 
ve _ WER Oe) 
O. = ty = nts = “90 x 103 = 0.47 rad 


Step 3: Prewarp the analog frequencies to yield the desired digital frequencies. 
wi, = tan (6/2) = tan(0.17) = 0.325 rad/s 
wi, = tan(@, /2) = tan(Ul27 0.726 rad/s 
Note: Prewarping is done to ensure the analog filter will ultimately 
yield a digital filter with the correct critical frequencies. This is best 
visualized in the following graph (Figure 2.1) of the relationship between 


the analog filter frequencies and the desired digital frequencies, 


eo — van oe 
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TABLE. 2.1 


ANALOG FILTER FREQUENCY TRANSFORMATIONS 


(after Ref. [1]) 


Filter Type or form 
to find 
ne, 
Lowpass filter 
1.0 
Uh p 
he, = — 
hate: 
0.0 w 
Highpass filter 
1.0 
om (he 
Bue crear 
0.0 w 
Bandpass filter 
1.0 
Wap ~ Wy 
hp, = 
Bue 
0.0 w 
Banastop filter 
1.0 
Bwas 
ape, = = 
~Wes ~ Wo 
0.0 w 
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$-form 

(to find H(s) trom Ayp,(s) 
replace $ in prototype 
with) 





Normalizing to the lowpass prototype, according to Table 2.1: 


w!, = 0.325 —+ w_ =1 rad 
Cc Cp 


wi, = 0.726 —> wa, = 0.726/0.325 = 2.234 rad 





Figure 2.1. Relationship Between Analog and 


Digital Filter Frequencies 


Nete: Normalizing is done to transform the prewarped critical frequencies 
of the analog filter to their relative equivalents in a lowpass prototype filter with 
a critical frequency of w,, = 1. This is done to enable the designer to take 
advantage of previously compiled tables or frequency response curves corresponding 
to the transfer functions of these lowpass prototypes. Thus, a prototype filter can 
be selected whose frequency response characteristics have the same shape as the 
filter that is being designed. The transfer function for the selected prototype is 
then converted to a transfer function that exhibits the desired frequency response 
characteristics using an appropriate substitution (Step 5). Figure 2.2, depicts the 


process of normalizing. 
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Bal Lowpass Prototype M Lowpass 





Figure 2.2. Nermalization 


Step 4: Find the order N of the lowpass Butterworth filter prototype. The order 


N to provide a gain of gp at an angular frequency w, is 


log 9 (10(—Mas/10) _ 1) 
2 logi9 Wa 
_ logy (10! — 1) (2.4) 
2log,9 2.234 
= 1.37 (Since N has to be an integer choose, N = 2) 


Ne= ‘Mas = —10dB 


Step 5: From Table 2.2 determine the transfer function for the lowpass filter. 














H p(s) = Hxrp,(s) = Hy p,(s) 
s= st ce ce 
(Fy 
eee (2.5) 
s?+V2s41 s= ss 
0.106 


3? +0.460s + 0.106 
Step 6: Determine the transfer function for the digital version of the filter using 


the bilinear transformation. 


tinee (3) 


OES (2.6) 
eo et | 0.068(= + 1)? 


ret me eeeode 0 c= dds + 0s13 
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TABLE 2.2 


BUTTERWORTH PROTOTYPE COEFFICIENTS 


aj 


1.0000 
1.4141 
2.0000 
2.6131 
3.2361 
3.8637 
4.4940 
5.1258 


(Table after to Ref. 4) 


a2 a3 

1.0000 

2.0000 1.0000 
3.4142 2.6131 
D.2o01 6.2361 
7.4641 9.1416 
10.0978 14.5918 
13.1371 21.8462 

Hz pp(s) = 


a4 


1.0000 
3.2361 
7.4641 
14.5918 
25.6884 
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as 


1.0000 
3.8637 
10.0978 
21.8462 


1 


a6 a7 ag 


1.0000 
4.4940 1.0000 
13.1371 5.1258 1.0000 


Tiras past eee 


Step 7: Obtain a plot of the digital filter frequency response to see if the design 


specifications have been met. 


Figure 2.3 shows that the design does indeed comply with the given filter 
specifications, in that the cutoff frequency is 0.628 = 0.2” rad and the stopband 


frequency is 1.257=0.4z7 rad with gains of —3 dB and —14 dB, respectively. 


Example 2.2 Chebyshev Filter 


Design a digital bandpass filter with the following specifications: 
e 1 dB ripple in the frequency band 600 to 900 Hz, 
e sampling frequency of 3 kHz, 


¢ maximum gain of —40 dB for 0 < f < 200 Hz. 


Step 1: Convert the critical analog design frequencies, w;, to the corresponding 


digital frequencies, 6;. 
Sampling frequency: f, = 3 kHz 


Lower ripple passband frequency: fe = 600 Hz => we = 12007 rad / s 
Upper ripple passband frequency: f= 00 Mz = 40, = 18007 rad /s 
Stopband frequency: f, = 200 Hz => w, = 4007z rad /s 


1200 
6 =e) = Wea war —e4t— 126 rad 
18007 
6 =w,T =w./f,= —C—ole d 
Sw Dy a) 3000 0.67 88 ra 


4007 
g! — ne — a 8 —— bs — < 
ew ny fe 3000 0.1337 = 1.54 rad 





Ripple band center frequency : 0% = \/0),6!, = 0.497 = 0.418 rad 


7X 
4 
4 


MAGNITUDI 
0.00.1 0.2 0.3 0.4 0.6 0.6 0.7 0.8 0.9 1.0 LA 1.2 





0.00000 0.62832 1.25664 1.88496 2.51328 3.14160 
FREQUENCG feels 


Figure 2.3. Frequency Response for Butterworth 


Lowpass Filter 
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Step 2: Prewarp the analog frequencies so that the desired digital frequency 


characteristics will be achieved. 


WU eatano,) 2) — tan 1.26/2) = 0.729 rad/s 
wi, = tan (0, /2) = tan(1.88/2) = 1.369 rad/s 
w. — tam. /2) —tan(0.418/2) = 0.212 rad/s 


Wy = 4/wyw!, = /(0.729)(1.369) = 1.00 rad/s 
Ripple band is: B' = wi, — w, = 1.369 — 0.729 = 0.64 rad/s 
From Table 2.1 convert the bandpass filter frequencies to the correspond- 


ing lowpass prototype filter frequencies. 


wip, =< ween we 
ee nae - 
W'BP WLPp 

w, 0.729 —1 

w, 1.369 1 

wi 0.212 —7.04 

w) 1.000 0 


Step 3: Determine the order of the Chebyshev prototype filter that meets the 
design requirements. 
For a lowpass Chebyshev filter the magnitude - squared characteristic is 


given by: 
: 1 
= Tr 2Ch lw) CR (w) (2.8) 


where C'y(w) is an N**- order Chebyshev polynomial and e indicates the degree 


Tolga (jw) 








of ripple. 
In this example the —40 dB gain translates to M = 10~? for the desired 


filter frequency response in the stopband, i.e., we want |H (jw)| <<) 10=- or 
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| H (jw)|” < 10~*. The design specification for 1 dB ripple in the ripple band 
corresponds to a value for € of 0.2589. 
Making an initial guess for the filter order of N = 2 yields: 








2 
il 
H rw Se a as sea EOS 
2(J ) il abe 0.2589 (2w2 — We We=7.04 (2.9) 
1 


= ————___—____,. = 4.01 x 107% 
1 + 0.2589 (99.12 — 1) 


which is not less than 107?. 


Proceeding further with N = 3 yields: 








3 
iL 
330) = == ee 
1 + 0.2589 (4w3 — 3wa)” |w,=7.04 (2.10) 
1 


= ——__________—_, = 2.044 x 107° 
1 + 0.2589 (1395.65 — 21.12) 


which is less than 10~*. Therefore, the design needs can be met using a third-order 

Chebyshev lowpass prototype filter. 

Step 4: Obtain the transfer function of the third order lowpass prototype filter 
with 1 dB ripple. Looking at Table 2.3, yields the following prototype 


transfer function: 


0.491 


a 2.14 
s3 + 0.9885? + 1.238s + 0.491 ( ) 


Hy p,(s) = 


Since this is an odd order filter, the constant 0.491 in the numerator was 
selected to make |H(j 0)| = 1. Recall that for Chebyshev filter design, the constant 
K in the numerator of the lowpass prototype filter is selected on the following basis: 

for N odd; 160) all 


(2.12) 


for N even; 260) = iaaie 
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TABLE 2.3 


CHEBYSHEV - PROTOTYPE DENOMINATOR POLYNOMIALS 


OWN 


wa 


(Table due to Reference 4) 
0.5-dB Ripple Chebvyshev Filter (e = 0.3493) 


5+ 2.863 

2 + 1.4255 + 1.516 

$3 + 1.25352 + 1.5353 + 0.716 = (5s + 0.62652 + 0.6263 + 1.142) 

$4 + 1.19753 + 1.71752 + 1.0255 + 0.379 = (52 + 0.3515 + 1.064) s2 + 0.845s + 0.356) 


g5 + 11725184 + 1.937452 + 1.309652 + 0.75255 + 0.1789 © 


(s ~ 0.3623)[(s + 0.1120)2 + 1.01162)[(s + 0.2931)? + 0.62522] 


$§ — J.159255 + 2.171854 -— 1.589853 + 1.171952 + 0.43245 + 0.0948 


((s + 0.0777)2 — 1.00852)[(s — 0.2121)? + 0.73822}[(s + 0.2898)? + 0.2702?) 


57 + 1.181256 + 2.412655 + 1.869454 + 1.647953 + 0.755652 + 0.28215 + 0.0447 


(s ~ 0.2562)[(s ~ 0.0570)? + 1.00642][(s + 0.1597)? + 0.8001 2}[(s + 0.2308)? + 
0.44792] 


s§ = |.]46157 + 2.656755 + 2.149253 + 2.184054 + 1.148653 + 0.573652 + 


0.1525s — 0.0237 
[(s + 0.0436)? + 1.00502)[(s ~ 0.1242)? + 0.85202)[(s + 0.1859)? ~ 0.56932] 
{((s + 0.2193)* + 0.19992} 


5s? + 1.142658 + 2.902757 + 2.42935 + 2.781555 + 1.611454 + 0.983652 + 


0.340852 + 0.09415 + 0.0112 
(s + 0.1984)[(5 + 0.0345)2 + 1.00402}[(s + 0.0992)? + 0.88292][(s + 0.1520)? + 
0.65532)}{(s + 0.1864)? + 0.34872] 


s19 — J 40159 + 3.149958 + 2.709757 + 3.440956 + 2.144255 + 1.527454 + 


0.627053 ~ 0.237352 + 0.04935 + 0.0059 
[(s — 0.0279)2 ~ 1.00332][(s + 0.0810)2 ~ 0.90512][(s -+ 0.1261)? + 0.71837] 
[(s -- 0.1589)* ~ 0.4612-)[(s + 0.1761)* +- 0.15892] 


rr Si esp 
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TABLE 2.3 (Continued) 


CHEBYSHEV - PROTOTYPE DENOMINATOR POLYNOMIALS 


mw he 


la 


1.0-dB Ripple Chebyshev Filter (¢ = 0.5089) 


s + 1.308 

s+ ~ 0.8045 ~ 0.637 

$3 — 0.73852 — 1.0225 — 0.527 = (s — 0.402) 52 ~ 0.3695 — 0.886) 

54 — 0.71653 — L.2S6s* — O.517s — 0.206 = (52 — 0.2105 — 0.928)(52 — 0.5065 — 0.271) 


53 — 0.706554 — 1.499553 — 0.693552 — 0.45935 ~— 0.0817 
(s — 0.2183)[(s — 0.0675) — 0.97382][(s — 0.1766)? — 0.6016] 

56 — 0.701255 — 1.745954 - 0.867053 — O.771Ss* — 0.21035 -- 0.0514 
[(s ~— 0.0370) -. 0.98172][(s - 0.1283)= - O.7187=][(s - 0.1783) — 0.2630] 

57 — 0.6979s® — 1.993555 = 1.039254 — L.lddds3 — O.382Ss- — 0.166615 — 0.0204 
(s ~ O.1853)[(s ~ 0.0346)* — 0.9867-][(s — 0.0968)= ~ 0.7912-j[ts — O.1399)- — 
0.4391] 

s§ — 0.696157 — 2.242356 — 1.211755 — 1.579654 + 0.898253 — O.35875* — 
0.07295 = 0.0129 
{[(s — 0.0265)* — 0.9898=]{(s - 0.0754)* — 0.8391 =][(s ~- O.1129)= — 0.56077] 

((s — 0.1332)- — 0.19692] 

59 — 0.694758 + 2.491357 — 1.383755 — 2.076755 — 0.856954 — 0.644553 + 
0.168445- -- 0.0Ss45 — 0.0051 
(s + O.1206)[(s - 0.0209)> -. 0.9919-][(s — 0.0003)2 — 0.8723*][(s — 0.0924)¢ 
+ 0.645742] [(5 -- O.1134)* — 0.34452] 

510 — 0.693759 = 2.740058 — 1.555757 — 2.626355 -— 1.188535 — 1.038954 + 
0.317853 — 0.1440s= — 0.02335 -- 0.0032 
((s — 0.0170)2 -~ 0.9938=][(s - 0.0767)= — O.7112"][(s + 0.0193)- — 0.8962-] 

[(s — 0.0967)* — 0.45672][(s — 0.1072)? ~— 0.15742] 





LA mw Wh ee 


2-dB Ripple Chebyshev Filter (e = 0.7648) 


s+ 1.965 

s? + 1.0985 — 1.103 

53 = 0.988s* — 1.2385 — 0.491 = (s + 0.49452 + 0.4905 — 0.994) 

5* + 0.95353 + 1.4545 + 0.7435 = 0.276 = (52 — 0.2795 — 0.987)(s2 ~ 0.6745 — 0.279) 


55 + 0.936854 + 1.68885) — 0.974452 + 0.58055 + 0.1272 
(s + 9.2895)[(s ~ 0.0895)2 + 0.9901 =][(s + 0.2342)2 — 0.61 192] 

5® + 0.928255 + 1.930854 — 1.202153 — 0.939352 — 0.30715 — 0.0689 
[(s — 0.0622)2 + 0.99342][(s ~ 0.1699)? ~ 0.72722][(s ~ 0.2321)2 — 0.2662 =] 

$7 = O.92315¢ — 2.176155 — 1.428854 — 1.357553 — 0.548652 — 0.21375 — 0.0307 
(s — 0.2054)[(s — 0.0457)2 — 0.9953=][(s — 0.1281)2 ~ 0.7982=}[(s - O.1851)2 — 
0.44292) 

s* — 0.919857 — 2.423056 + 1.655253 + 1.836954 — 0.846853 — 0.457852 — 
0.10735 — 0.0172 
[(s — 0.0380)2 — 0.99652][(5 — 0.0997)2 — 0.8448=][(s + 0.1492)? — 0.56442] 
[(s ~ 0.1759)? — 0.19822] 

5? — O.9175s* + 2.670957 — 1.881556 — 2.378155 — 1.201654 — 0.786353 + 
0.24425* — 0.0706s — 0.0077 
(s ~ 0.1593)[(s ~ 0.0277)2 — 0.99722}[(s + 0.0797)? + 0.87692][(s — 0.1221)? + 
0.65092)}f(s — 0.1497)? — 0.34632] 

510 — 0.915959 — 2.919858 ~ 2.107957 = 2.981556 — 1.683053 — 1.243554 + 
0.455453 — 0.182532 — 0.03455 — 0.0043 
[(s — 0.0224)2 — 0.99782}[(s ~ 0.1013)2 — 0.7143=}{(s ~ 0.0651)= — 0.90012) 
[(s — 0.1277)2 — 0.45862}[(s — 0.1415)2 — 0.1580] 


EO a ete ee eeenpenasee> 


Step 5: Convert the lowpass filter prototype to a bandpass filter prototype. 
From Table 2.1: 


Hpp(s) = Axi pp(s) 


a2z4w2 
s=—,z,*% 
0.491 
= ae See el (2.13) 
(aa) 1F 0. 988 ae ah)? a 1.238 (ot ~+1 ) oF 0.491 
0.129s° 


~ 38 + 0.63255 + 3.507s! + 1.3945? + 3.5075? + 0.6325 + 1 
Step 6: Determine the digital version of the transfer function using the bilinear 


transformation. 


Hpp(z) = Hpp(s) 
s= eat 


7 0.19 (#51) 


ar cae sn tsa) 4 3.507 (1) 


= 








+1238 (2 <=>) +0491 
+] 
_ _ 0.011 (2° ~ 3z* + 32? — 1) 
28 + 2.15324 + 1.7862? + 0.545 
(2.14) 
Step 7: Confirm that the filter design meets specifications by obtaining a fre- 


quency response plot. Figure 2.4 indicates the specifications for a pass- 


band of 1.26 to 1.88 rad with 1 dB ripple have been met. 


Example 2.3 Elliptic Filter 
A lowpass elliptic filter with the following specifications is desired: 
e passband ripple of 0.5 dB, 
e passband ripple-edge frequency of 2 kHz, 


e stopband gain should be at most —20 dB for frequencies greater than 4 kHz, 
and 


e sampling frequency is 20 kHz. 
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3.14160 


N= 6 


Step 1: Convert the critical analog design frequencies to digital frequencies. 


Step 2: 


Step 3: 


Step 4: 


Sampling frequency: jo AU ele 
Passband ripple-edge frequency: f; =2kHz => w; =4 x 10°m7 rad /s 


Stopband ripple-edge frequency: f2 =4 kHz => w2 = 8 x 10° rad /s 


4x 103) a 

_— => _ (4x 10°) = 

ee se 0.27 = 0.628 rad 
(8 x 10°) x 

Gy = a = 02] fs gas = 0-8" = 1.88 rad 


Prewarp the analog frequencies in order to determine the appropriate 


lowpass prototype filter. 
evan 2) = Udo tad/s 


= tamjo,/2) = 1.376 rad/s 


Determine the order of the elliptic prototype filter that meets the design 
specifications. Find the ratio, R, of the stopband frequency, w5, and the 
passband frequency, w}. 


1.376 


R= wil = 7355 


= 4.234 (2.15) 


From Table 2.4b determine the prototype filter order and the resulting 
transfer function. Looking at Table 2.4, it can be seen that a filter of 
order N = 2 has a value for R of 2.76261, which is more than sufficient 
to meet the design specifications. 

Obtain the transfer function of the second-order lowpass prototype filter. 
Again, looking at Table 2.4, the prototype transfer function is determined 


to be: 


0.1s* +- 0.5338 


H SS Sea a ae 
LPe(S) = 39 8094s + 0.5607 


(2.16) 
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TABLE 2.4a 
ELLIPTIC FILTERS 
Generalized Transfer Functions 


(due to reference 5) 





Nie 
s? + Ag) Hos? + Hy Aoi 
H2(s) = Hp) —-——__—_- = = 
2(s) (22 eon s*+Bys+ Bo 
N=3 
Ho s* + Ao 
H;(s) = —— -_ —_——_——_——__————__ 
a(s) stso s*+Biys+ Bo 
Hos* + Ho Ao 
H3|3) =) ee Oe ee 
s3 + (By + 89) 8? + (Boi + Biiso) s + Boi $0 
N=4 
i). -_ “ao eee 
: ” 8? + Bus Boy 82 + Bios Bap 
Oe Hos* + Ho (Aoi + Ao2) 8? + Ho Aoi Ao2 
s4+ (Bi, + Biz) s? + (Boo + Bir Bi2 + Bor) s? + (Bi1 Bor + Boi Biz) s 
+ Bo: Boo 
(Nes) 
Reto Hys* + Ho (Aoi + Ao2) 8? Ho Aoi Ao2 


s° + (By, + Byo +39) s4 + [Boo + Bi: Bi2 + Bor + Biiso + Bizso] s3 
+ [B11 Boo + Bor Bi2 + Borso + Bi Bi2s0 + Boi so] s? 


+ [Bo1 Bor + Bi Bo2s0 + Boi Biz $0] s + Boi Bor so 
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TABLE 2.4b 


ELLIPTIC FILTERS 


t 


(due to reference 5) 


Ta 


Dvs) “| 


+5.  Nodd =| 


Veven 


(a) Passband rippie 0.5 dB, stopdand gains. - 20, -30, -40, -50. -60. -70 dB 








Passband mpple = 


05 dB. slopband gain = -20d8 





N t Ay By 8,, H4/'So it, 
2 1 S$ 33789 0 S6hoH0 0 809390 0 100220E + 000 2.76261 
3 | | 7S640 0 808321 0 359160 0.2062) 4F + 000 1.42189 
0.667292 

i" =) 4 34105 0 G11198 1) 931989 0.}00219E + 000 1.43188 

a 1 21841 0) 927332 0.136543 
3 | | 65076 0) 827787 0 412816 0. JN389SE + 000 1.04465 

EY 107211 0 973640 0 1939S 0 667292 

\ 4 36790 0 611899 U 933855 0.1002} 8E + 000 1.01853 
6 2 1 19243 0 933830 0.1$6224 

a | 02486 0.990620 1.017576 

\ 1 YIN 0 828092 0 413652 0) 303861E + 000 1.00545 
7 2 1 06401 () 976379 1) OS6384 0 667292 

3 | 00870 0 9966x | 1) UU6219 

l 4 36811 0 411836 933868 0. LOOL92E + 000 1.00192 
2 1 19207 () 933928 156548 

A Oss 0 991634 020051 

a | 00306 0 998827 0.002197 

1 | 64927 1). 828057 0 413695 0.303786E + 000 1.00068 
& 2 1 06190 0976512 0 056505 0 467292 

1 | (0778 O 997051 0 ON7093 

s | 00108 0 9995K6 {} WOO77§ 

Passband ripple = 05 dB; stopband gain = -30 d8 
N ; A, 8, 8., Ho/So 2, 
2 1 9 $1248 0.318702 0 639007 0.316294E - 001 4 80880 
3 | 2 46997 0 597384 0.382044 O.I2E878E + 000 1.92322 
0 $1176 

‘ | 6 46603 0 398996 7822201 0 316297E - 001 1.32446 

2 1 47114 0.798764 0 191032 
A 2 14490 0 648724 () 480774 OV M8O7E + 000 1.42942 

a 1 18132 0.907216 0 088080 U SI1761 

\ 6 38228 0 402050 O 828822 0.316296E - 001 1.05394 
6 u 1.38680 0 82682} 0 237025 

4 1 074374 0.988727 0 039181 

I 2 13439 0.650591 0 483325 0.11 8701E +000 1.02299 
z 2 115171 0 920785 0 108319 O $1176 

3 | O3168 0 981925 0.017159 

! 6 37941 0 302154 O 829082 0 3162 89E — 00] | 00989 
: 2 1 38394 0.827819 () 23%663 
‘ 1 1 06d01 O 9648 9x 0 048043 

4 | 01359 0 992137 0 007464 

I 2 13409 0 6506.36 ) 484554 O 1 E8683E + 000 1.004827 
9 2 | 1$070 0 921256 0.109143 0.511761 

3 1.02680 0.984652 0.021005 

4 | OOS86 0.996589 0 003238 
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TABLE 2.4c 
ELLIPTIC FILTERS 


(due to reference 5) 





Passband rippie 1 dB, stopband gains - 20. - 30, - 40, -50. -60. - 70 dB 


Passband ripple 


= 1 dB: stopband gain = - 20 0B 

















N t AN By, (eh, Ho: So Q, 
? 1 fy SRR 0 a97234 () 676727 O. LOOT SSE + 000 2.32474 
: | 1 5856S 0.790229 0) 242927 0.28 LOXOE + 000 | 30797 
O S6S168 

; | A wKI47S 0.536633 0 768217 0 LOUL8SE + 000 1.09029 

2 1 15956 0.926578 0.099029 
. | 1 SIKS2 0 #08039 0.318242 Q 279829E + 000 | 02426 

a 1 OdkK6 0 978703 0 032771 0.565 16x 

| 3 80X73 0.537071 0 769217 0 1O0L84E + 000 1.00902 
) 2 1 14367 0 YAI19S 0.110761 

3 1 O1SS0 0.992101 0.010684 

! 1 SI782 0.808242 0.31627 O0.279R SE + 000 1.00290 
? a | O32] 0 977941 0 036573 0.565 16% 

3 1 00397 0 997446 0.003444 

1 + 80K66 0.537076 0 769228 0.10018SE + 000 | 00093 
i 2 1 14380 0 933266 O. | 10888 

| 1 O1405 0 992814 UO OLLSsl 

4 } 00160 0 999176 0.001111 

} 1 SIK12 Q RO&OOY 0 318714 0.27961 1E + 000 1 00030 
2 2 ) $421 0 9779 %6 0 036644 0 S65168 

3 1 O88) 0 997680 GO OUIKIS 

4 1 uoos? 0 999734 0 0OU3S9 

Passband ripple = 1 dB, stopband gain = - 30 dB 
N i Ay By 8,, Hoi So (Yh 
Ry | 7 KB1SX 0.279699 0 $38632 0 316284E - 001 400423 
; | 2 20293 0 546596 O.ALI98I O 133223F + 000 } 73254 
Q 430700 

5 l Sir2or2 O 348111 0 682880 0.316286F - 001 1.25040 

2 | 37628 0 03477 0 148347 
A | 1 96914 0 633529 0.38% 3990 0 LELTSOE + 000 1 G9SS84 

2 | 391K 0 9144045 0.064612 0 4330706 

1 § 67741 0 350184 0 687078 0 316278E-001 1 03799 
6 2 1 31626 0) ®2%467 0.179733 

‘ | OS46% 0 930K 0027111 

] 1 96322 U 6 148R4 0 VROOSS QO LITL236. + 000 1.01536 

= 1 11870 0 925872 0 077672 0) 430700 

3 1 02200 0 9%S162 0 Q11201 

} 5 67614 0.350237 0 587189 0.31626SF - OUI 1 00625 
‘ 2 1 31463 0 429171 0.180623 

3 | 04690 0 969003 0.032478 

4 | 00893 0 993908 0.003594 

| 1 96311 0 644908 Q AROL17 0 1ILIO9E + 000 | 00255 
9 2 | 1k 0 926187 O O78O83 0 330700 

3 1 OV89T 0 987211 0 013399 

4 1 00168 0 997805 Q OOLBRY 
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Pasaband nppie 2 dB, stopoand gains: ~ 20, -30, ~ 40. ~50, -60, ~70 08 


TABLE 2.4d 


ELLIPTIC FILTERS 


(due to reference 5) 








PassbDand ripple = 2 dB: stopband gain = ~ 20 dB 
N i As Eh 8,, Ho/So 2, 
2 | + 8094 | 0 4649] 0) $17326 0 LO0L03E +000 1.94332 
; | 1 42919 0 793180 0) 2048089 0 284443 E +000 1 20808 
() 458898 

3 | 3 25882 0) 466218 0) $9726 VU. 100 1L02E * 000 1.05569 

2 | 10768 0 Y3SS6s 0 06348§ 
‘ I 1 39114 QO #07314 1) 223995 0 283878E + 000 | 01867 
ss 4 102976 0.9% 1070 1) OF 8680 () 458898 

! ¥ 28087 () 4864337 0 §97479 0, LOQLOZE + O00 1 0457 
6 2 1 My ys 0 940286 0) 069417 

i 1 O0¥45 () 994532 1) OOS 394 

] 1 (W096 0) XO7 M4 ) 224146 0 2518 YE + O00 1 OOL28 
3 y 127490 () OM 244 0 0203462 () 45% 898 

i j 242 () Y98428 ) OO1SSt 

| 3 257346 0.444314 0 597662 1), LOOOSNE + 000 1.00037 
" 2 19918 0 940278 0 069488 

3 1 00782 0.994918 1). 005883 

4 | 00069 0 999548 0 000846 

| 1 19187 () $0705 0 224311 Y ISAIWE + OO) | OOOH 

S 1 ay7ss (2962537 0 020319 1) 3SKMUK 
uf 

ty | tmo224 1) GOX S40) 1) LO) 694 

a 1 00020 Ob QY9X70 1) OOOL29 

Passband nppie = 2 dB. stoppand gain = -30dB 
N ' A, Bn 8,, H9/So i. 
2 6 42917 Me TSSOTS () 433983 O.3S162S59E - OO] R29 235 
1 1 482%) (0) $Q377% 1) 238474 1) JOM S3E = 000) 1 SS040) 
() J4S92% 

i | 4 9944K 1} V2) 84 ) $7291 U 31A2S8F - 101 1 18280 

uy | SaaS 0 X200246 0) 105666 
= | 1 798% 0 634366 0 285631 0 JOI947E + 000 1 06594 

S] | e838 0 927270 0 032702 () 345928 

| 4 967AS Ow 1 S%9SA2 0 316259E = O01 1 02460 
6 eS Lh cleat 0 840500 0 124608 

‘ 1oksrsy 0971673 0) 016632 

| 1 *9218 Q) 638340 0) 286656 0 102920F +000 1 00929 
os Bs 1 Ost8] () 93S9RKS 0 980009 1) 345928 

i to1dol 0 OX¥EIS 1) OO6386 

i 496715 0.313454 0) §3960% O VWa28 hE - 001 1 0383 
: = | Janus 0 ¥409%3 0 125007 

: | VS4 () 97S] 37 1) 019424 

4 J} mS) 0) Y9SH47 0 002438 

| | “9236 i Aisr1o () 280689 0 1028900 «000 1 iotas 

3 | Ux736 0 961S0 9 OS017% i) $4592 

3 it dy) eed 0 OX WAG 0 007482 

a Poso2 () YOK417 0 OOUG AD 
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Step 5: Find the analog transfer function of the desired filter based on the lowpass 
prototype. 
This is done by finding a scaling factor, a, to convert the lowpass proto- 
type’s cutoff frequency, w;,, to the desired cutoff frequency, w}. 


For the prototype: 
R=2.76261 ; wip =1/VR =0.6016 


For the desired filter: 


Therefore, the scaling factor, a, is: 
a=w,,p/w, = —— = 1.85 (2.10) 


Thus, the actual transfer function, Hz, p(s), can be found as follows: 


| 


Hy p(s) = Hzpp(s) - = Hyp,(s) 


0.1s? + 0.5338 
s? + 0.80948 + 0.5667 | ,_, 35, (2.19) 
0.1(1.85s)* + 0.5338 
(1.85s)? + 0.8094(1.85s) + 0.5667 
_ 0.34238? + 0.5338 
3.4238? + 1.497s + 0.5667 


Step 6: Determine the digital version of the transfer function using bilinear trans- 








s=1.853 





formation. 


H, p(z) = Hy p(s) 
= Se 
_ 0.3493(2 =1/2-F 1) discas 
3.423 (z—1/z+1)° + 1.497(z — 1/z +1) + 0.5667 (2.20) 
_ 0.87612z? + 0.383z + 0.8761 
~ §.48722 — 5.7126z + 2.493 
_ 0.1597 z? + 0.0698z + 0.1597 
~  -.z2 — 1.04112 + 0.4543 
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Step 7: Verify that the design meets the specifications from a frequency response 
plot (Figure 2.5). 
2. Summary 
This concludes the review of the use of the bilinear transformation to 
design recursive filters based on Butterworth, Chebyshev or elliptic analog pro- 
totypes. As can be seen, this method is very involved in terms of algebraic ma- 
nipulation. The following sections will introduce the direct design technique. It 
will be shown that the direct design procedure reduces considerably the amount of 
algebraic calculations required, and eliminates the somewhat confusing procedure 


of prewarping. 


C. DIRECT DESIGN TECHNIQUE 

The direct design technique is based on the bilinear transformation in the 
following manner: 

The tables of coefficients for analog lowpass prototype filters, i.e., Butterworth - 
Table 2.2, Chebyshev - Table 2.3, and elliptic -Table 2.4, were converted to 2- 
domain versions through the use of the bilinear transformation, and can be shown 
to have a critical frequency of 7/2. The reason, as stated in Section A, is that 
the analog prototype filters comprising these tables all have a critical frequency of 
w. = 1. For the bilinear transformation, the following substitution for s is used in 


the analog prototype transfer function to convert it to a digital prototype transfer 








function: 
—1 1 
ae ones = (221) 
Zan Il 1l—s 
An analog critical frequency of s = jw, = 1 implies: 
ee) 1+) 
ee (22D) 


7 1 40 7 l= 
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since z = e/°, this gives: 


Sipe a = 1I7/2 (2993) 
ts 
Thus, it can be seen that an analog critical frequency of w, = 1 yields a 


digital critical frequency of 9, = 7/2. Tables of prototype filters were thus created, 
specifically, Butterworth - Table 2.7, Chebyshev - Table 2.8 and elliptic - Table 2.9. 
The use of these tables to design a digital filter based on analog specifications 
involves the following generalized steps. To further illustrate the technique, a 
detailed description of its derivation and use, including examples for the various 
filter types, will be presented. 
1. Direct Design 
Step 1: Determine the appropriate desired filter type based on the design speci- 
fications, i.e., Butterworth, Chebyshev or elliptic. 
Step 2: Convert the analog filter frequency specifications to digital equivalents. 
Step 3: Use the digital-digital frequency transformations of Table 2.5 to normal- 
ize the desired filter’s design frequencies, so that the appropriate lowpass 
prototype filter may be selected. 
Step 4: From the normalized design frequencies obtained in Step 3, determine 
the lowpass prototype that meets or exceeds these requirements. 
For Butterworth and Chebyshev filters - this is done using the design 
curves of Figures 2.8 and Figures 2.10 - 2.12. 
For elliptic filters - Table 2.9 is used. 
Step 5: Obtain the actual filter transfer function from the lowpass prototype, 


based on the digital-digital frequency transformations of Table 2.5. 
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TABLE 2.5 
DIGITAL FILTER FREQUENCY TRANSFORMATIONS 


(Due to Reference 2) 


Type Transformation Design Constants 
(Replace z in LP digital 
prototype with) 


_ sin( 6, /2-8, /2) 

















Lowpass == OS 
Pp isaz sin( 6, /2+6, /2) 
cos( 8 2-6, 2) 
2. Highpass == ee a ’ 
s( ¢ [248 /2) 
2_ 2ak k-1 f / 
_ 6 /2+6%, /2 
: z ot a cos{ Y,, c4 
3. Bandpass oe eer ee ae FASE ae (61, /2—6% /2 
eae eee AGETET) 
0 
= iol ! = 
k = tan Scot (@,/2 — 6,/2) 
2 2a 1-—k d / 
= as cos( 4 /2+8 /2) 
4. Bandsto <— a ; 
P 1 Pie e & co5(8%,/2—8, /2) 


k = tan fo tan (0, /2 — 0/2) 
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As can be seen from the previous design summary, the crux of this method 
is the use of the digital-digital frequency transformations of Table 2.5, which merit 
explanation [2]. 

These transformations enable the user to convert the lowpass digital pro- 
totype, Hz p,(z) (be it Butterworth, Chebyshev or elliptic) to the actual lowpass 
(LP), highpass (HP), bandpass (BP), or bandstop (BS) filter transfer function. 

The transformations provide a means of transferring the stability, inherent 
in the prototype filter, to the actual filter; that is, the poles of the actual filter will 
lie inside the unit circle, as they do for the prototype filter. For this reason the 
frequency response of the lowpass prototype filter at a specific frequency of 9, must 
be the same as the desired value of the frequency response of the actual filter at 


its corresponding frequency of 6/. 
Hips (€7%) = Ar ype (ci*) (2.24) 


where, 
6, = prototype frequency 


6 = desired filter frequency 
iijije— Ul. np. DP aor BS 
For example, when transforming a lowpass prototype filter to a desired 
lowpass filter, one wishes to maintain the integrity of the magnitude characteristic 
of the prototype filter, while expanding or compressing the frequency scale so that 
the critical frequency is changed from the prototype value of 6, = 7/2 to the 
critical frequency of the actual filter #,, (see Figure 2.6). 
According to Table 2.5, this involves the following substitution for z in 
the lowpass prototype transfer function. 


a (2.25) 


Py Sie 
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M Lowpass Fretotype Lowpass 





Figure 2.6. Lowpass Prototype - Lowpass Transformation 


For the critical frequency of the prototype, 9,, to map to the critical 
frequency of the actual filter &, the following must be true: 

Substituting, e?%c = z on the left in Equation (2.25), and e/%> =z onthe 
right gives, 


- of 
ec aay 


f 
EL 


Ps 
1 — ae 


(2.26) 


and solving for the design constant, a, yields, 
eo sin (6, /2 — 61/2) 
sin (0, /2 + 61/2) 
Similar arguments apply to HP, BP, and BS filters. 
Thus, the transformations of Table 2.5 accomplish the appropnate map- 


ping from the lowpass prototype digital filter to the actual filter, and ensure that 


stabilty is maintained. 
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2. Butterworth Filters - A Direct Design Approach 


The design curves shown in Figures 2.8a - d, describe lowpass prototype 
filters with a critical frequency of 6, = 7/2, and of order N = 1 through N = 8. 
Higher order filters can be realized by cascading the appropriate lower order filters 
characterized by these curves. 

The curves were constructed by using the s-domain version of the But- 
terworth filter coefficients for N = 1 through N = 8, and then determining the 
z-domain version of these transfer functions through the use of the bilinear trans- 
formation and the conversion tables obtained from Reference 3, (see Table 2.6). 
These tables establish the necessary relationships between the coefficients of an 
analog filter and its digital counterpart. 


a. Generalized Analog Transfer Function 


= Ao Jo Ale g te sige +... 4 Aps* 


RS a 2 
(s) Bo + Bis + Bos? +...+ Bys* ot) 
b. Generalized Digital Transfer Function 
ay =) —k 
H(z) = COMO aegis Ses (2.28) 


1+ 6:27! + bgz-2 +...4+ bpz-* 


For example, in the case of N = 1 the bilinear transformation digital 


filter coefficients, in terms of the analog filter coefficients are: 


A Bo + Bi 

ao (Ao + A,)/A 
ai (Ao = A,)/A 
b (Bo -B,)/A 


Note: C is the critical frequency, w, = 1. 


For a first-order Butterworth prototype filter described by: 


2a) = — (2.29) 


4] 


TABLE 2.6 


BILINEAR TRANSFORMATION DIGITAL FILTER COEFFICIENTS 
IN TERMS OF ANALOG FILTER COEFFICIENTS 


A 

6 
A B+ B,C ° 

Oy 

02 
0, | (Ag-4yC)/4 

by 
by | (Bg-84C)74 

b2 


(Ist order) 


& | Bo+B)C+92C2 +8505+84C% 

Og | (Agra. CedgC2eayCreagcty/A 
0, | (449+24,C-245C5-6440%)/8 
02 | (649-2462 +64 407) 74 

Oy | (449-24, C +2409 -344C41 70 
Og | (Ag-h,Crhgc? -dyCoagc* 7A 
b, | (485+28,C-285C2-48,C7! 74 
b> | (689-2820%+684C7 1/4 

by | (485-28,C+283C°~s84C7) 74 
(Bo-8, C+B,C*-34C7+8gC41/4 


oO 
a 


(4th order) 


the values for the aigitai filter 


The transfer func 


H(z) = 


(due to reference 3) 
A |185+8,C+8,C7+84C> 


Bo+8,C +8267 Op | (AgrAyCraaC®+A4C5)/4 
(Agr AyCoA.C2)/4 ~ | (Sages C-agC2— 30503) 74 
(2A9=249C°)/A 02 | (3Ag-0,C-a9C2+3A5C7) 4 
(Ag -AyC+anc2)/a 0y | (Ag-AyC+agC2=A3C7) /A 
(2869-28207) /4 by | (389*8,C-82C%-38.C7 1/4 
(Bg-8,C #820714 by | (389-8,C-82C7+385C7)/48 
by | (Bg-8,C+B2C%-8xC7)/4 


(2nd order). 
(3rd order). 


A | B48, C+82C7+8yC5+84C*+85C° 
(Agra, Cob c2eAsC3 +A gC* oOgC>)/4 
oy | (Sag +34. 64h C2 -04C>- 2h go*-Sa5C°) JA 
(OA +24, C-2d9 C220 5C> +20 4C* + 100507) 74 
(1049-2, C-2A 907 24403420404 = 104507178 
04 | (SAg- 3A, C#haC2 #4 4C3- 344644 SAQCPI/A 
(Ag, Cage? =0503+44C*-A507) JA 
by | (S85*38,C+82C2-8,C3-38404-55gC7178 
b> | (1085 +28,C-28 2627-28503 +284C*+1085C7)/8 
bs | (1C89-28,C-282067 +28 40> +2RgC*~ 108407) /4 
bg | (89-38, C+BC° +85C°-384C*+585C7)/A 
be | (Bg-8,C+82C*-85C7+84C*-B5C7 1/8 

(Sth order) 


coeficients are: 


A 1+1=2 
a, (1—0)/2=0.5 


tion H(z) is therefore: 


0.54+0.527! 0.52 + 0.5 


1 => Ay (z) = way t (2.30) 


A plot of this transfer function is illustrated by the curve (N = 1) in Figure 2.7. 


42 


3 
v) 


MAGNITUDE 


0.6 0.7 0.8 0.9 1.0 1t.1 


U.5 


0.2 0.3 U4 





0.00000 0.62832 1.25664 1.88496 2.51328 3.14160 
FREQUENCY, THETA 


Figure 2.7. Frequency Response for Butterworth 


Lowpass Filter 


43 


Table 2.7 is a summary of the analog and digital versions of lowpass prototype 
Butterworth filter transfer functions for filter orders up to N = 8. 
To illustrate the direct design approach the design problem of the 
previous section will be used (Example 2.1), where the problem statement was: 
Design a digital filter for a 20 kHz sampling rate that is flat in the 
passband of 0 to the 3 dB cutoff frequency of 2 kHz, and has a gain of not more 
than —10 dB for frequencies greater than 4 kHz. 
Step 1: A Butterworth design is called for because a flat passband is desired. 


Step 2: Convert the critical analog design frequencies to digital. 


Sampling frequency: f, = 20 kHz 
Cutoff frequency: fo =2kH2 == wo, =4x 107 


Stopband frequency: f, =4kHz => wy, =8 x 10?m rad/s 


- 7 (4 x 10°) x 

oe = 1) ee = Wyma = 0-27 rad 
(8 x 10°) x 

ip = Pd ie — Wee = Doon = 0.47 rad 


Step 3: Since the “design chart” is based on a critical frequency of 9, = 7/2, 
digitized design specifications need to be normalized to this frequency, 
so that the order of the required lowpass prototype filter can be deter- 
mined. For this process the following frequency transformations apply 


(from Table 2.5, prime denotes desired filter). 


e Prototype filter from desired filter: 


at 
on ec) — a 


~ l= aei% 


(2.31) 
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TABLE 2.7 


LOWPASS PROTOTYPE BUTTERWORTH FILTER 


Filter Order 
N 


oe -®- WW b> 


TRANSFER FUNCTIONS 


(due to reference 4) 


Ht p,(s) 


1 
s+1 


| 
37+1.41423s4+1 





ee 
33+2324234+1 


1 
34#4+2.613139+3.414237+2.61313s+4+1 


1 
35 +3 .236184+5.236139 


25 .236i65 +3 .2261s+1 


1 
3° +3 .863735+7.464134+9.141633+7.4641 32 


+3.86373+1 


1 
3’ +4.49405°+10.0978s°+14.59183% 


+14.59183°+10.09783s7+4.4940s+1 


i 
3®+5.125887+13.13173%+21.846285 


+25 .6884 37+21.846238°+13.13713s" 


+5.12583s+1 
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Jes) 


0.5(2+1) 


z 


z+1 * 
3.41224+0.59 


(z+1)° 
623422 


(z+1)? 
10.6383244+5.17 22 


z+1)° 
18.9425+11.99623+1.0549z 


z+1 8 
33.797 25 +26.28424+4+3 .86227+0.0592 


: (z+1)" 
60 .367z7+55.5372z5+11.633 z9+0.4624z 


z+1)5 
107.896z°4114.4382%+31.496z4 


+2.162z7+0.008 
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3.14160 


e Desired filter from prototype filter: 


Bld 
95 oe) ee 


eo = Taree 
where (232) 
_ sin (6,/2 — 4 /2) 
sin (6,/2+ 4, /2) 


g'. and @, are the critical and stopband frequencies (respectively) used to 
determine the optimum prototype filter order. 


Equation (2.31) is applicable to this problem, and thus: 


a 

0 = 0.40 
Ce = Ring 
G, = 0.20 


_ sin(0.67/2 — 0.27/2) 
~~ sin(0.57/2 + 0.27 /2) 
_ sin(0.157) _ 0.454 


= = (Ea 
sin(0.357) 0.891 


Therefore, 
ah = eJ9-4™ _ 0.510 
~ 1—0.519e70.4" 


0.97¢21-779 5 atG 
Vo7enaaT (2d) 


= 6 = 2.349 rade— ion enad. 


Step 4: Once the normalized stopband frequency is obtained, (in this case 0.757), 
the design chart is used to determine the lowest order filter that will give: 


— 3dB at 0, = 0.57 


less than — 10dB at 0, = 0.77% 
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Turning to Figure 2.8, it can be seen that the curve for N = 2 meets 
these requirements. At this point, it is worth noting that a second-order Butter- 
worth filter was also determined to be required using the traditional approach in 
the previous section. 

In Table 2.7 it can be seen the prototype low-pass filter transfer func- 
fon tor N = 2 1s: 


(Ze) <i roan Sy S| 


per) = 3.4122 + 0.59 = 3.4122 + 0.59 


(2.34) 


Step 5: Once the prototype filter is obtained from the design curves, it is neces- 
sary to determine the actual transfer function that will meet the design 
specification, such that 6, is 0.27 and @, is 0.47. 

Using the digital filter frequency transformations of Table 2.5, the 


transfer function is determined as follows: 


Hy p(z) = Hripp(z) 














— ieee 
_ ge 250 ck I 
~ 3.4122 + 0.59 z= Ea 0.510 
27 22 + | 


= 2.35 
14.8422 — 17.02 + 6.14 ( ) 


which is the same transfer function obtained in the previous section. 

Step 6: A computer generated frequency response plot of the filter is used to 
verify that the design fulfills the design specifications; Figure 2.9 reveals 
that this is indeed the case. 

D. LOWPASS PROTOTYPE TO HIGHPASS, BANDPASS BAND- 

STOP FILTERS 
Once the lowpass prototype filter is obtained, a highpass, bandpass or band- 
stop filter can be derived using the digital frequency transformations of Table 2.5. 


In the next section, an example is presented that illustrates this process. 
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1. Chebyshev Filters - A Direct Design Approach 


As stated in the introduction to this section, this design approach is of 
course applicable to Chebyshev filter design problems. Again using the bilinear 
transformation techniques cited in [2], the digital filter coefficients corresponding 
to the analog filter coefficients for prototype Chebyshev filters were determined. 
The frequency responses of the resulting transfer functions listed in Table 2.8 were 
then plotted to obtain a set of design curves with a cutoff frequency of 7/2 (see 
Figures 2.10 thru 2.12). 

The design problem of the previous section will be used to illustrate ap- 
plication of this new technique (Example 2.2). 

The problem statement called for a Chebyshev digital bandpass filter to 
be designed to meet the following specification: 

e 1 dB ripple in the range of 600 to 900 Hz 
e Sampling frequency is 3 kHz 
e Maximum gain of —40dB for 0 < f < 200 Hz 
Step 1: Convert the critical analog design frequencies to digital: 
Sampling frequency: Jee Shale 
Lower Passband frequency: fe = 600 Hz 


Upper Passband frequency: f, = 900 Hz 


Stopband frequency: j= 200nz 
6; PS 7) fy ee = Ode = 1:26 rad 
ae Vee We} ls 3000 = U. = 1. 
2(900)zx 
g! = ved ky =e) = |. 
=) walt 3000 0:67 = 1-88 rad 
2(200)x 


Cae ee — = loot — Ole pad 


3000 
Ripple band center frequency: 0) = ,/0,6, = V2.35 = 1.54 rad 


o3 


TABLE 2.8 


LOWPASS PROTOTYPE CHEBYSHEV FILTER TRANSFER FUNCTIONS 


(due to reference 4) 


1/2 - dB ripple (ce = 0.3493) 


H(s) 
2.863 
3+2.863 
1.43 
3?+1.4253s+ 1.516 
0.716 
934+1.253827+1.5358+0.716 
0.944 





34+1.19733941.717324+1.0255+0.379 


H(z) 


—— 


2.863z+2.863 
3.863z+1.863 
1.43z74+2.86z+1.43 
3.941z2+1.032z+1.091 
0.7162°+2.148z7+2.148z+0.716 
4.504z9—0.570z2+2.360z—0.566 
0.358241 .43229+2.1482°+1.432z+0.358 
5.318z4—2.828z9+4.840z2—2.140z+0.874 


1 - dB ripple (e = 0.5088) 


a0) 


(s) 
1.965 
3+1.965 
0.983 
37+1.098s+1.103 
0.491 
93+0.9883s2 +1.23834+0.491 
0.891 


Pon e |S 


al 


947+0.9539341.4543924+0.7438+0.276 


H(z) 


1.965z+1.965 
2.965z+0.965 


0.983z7+1.966z+0.983 
3.201z27+0.206z+1.005 


0.49129+1.4732741.473z4+0.491 
3.717 g2=1 077 ee ea oe 


0.246 2*+0.984z°+1.476z7+0.984z+0.246 
4.42624—3.31629+4.748z2—2.476z+1.034 


2 - dB ripple (€ = 0.7648) 


H(s) 
1.308 
s+1.308 
0.506 
397+0.8043+0.637 
0.327 
33 +0.73832+1.02238+0.327 
0.164 


Bone fe 


34+0.71633+41 .25632 +0.5173+0.206 


o4 


7) 


1.308z+1.308 
2.308z+0.308 
0.506z7+1.012z+0.506 
244122 —0.726240.833 
0.32727 +40.981z7+4+0.981z+0.327 
3.087 29 —1 1352040 O21 0957 


0.164 24+0.656z°+0.984274+0.6562+0.164 
3.69524 =3 S742944.79422—7 778s) 5 


6 


0.045 
7 8° 41.151894+2.4138° +1.8698* +1.6488°+0.7568* +0.2828+-0.046 


TABLE 2.8 (Continued) 


LOWPASS PROTOTYPE CHEBYSHEV FILTER TRANSFER FUNCTIONS 


1/2 - dB ripple (¢ = 0.3493) 


H(s 


79 
09 1.17384 1.937 8541.31082+0.753840.179 


9.09 


7°4+1.15504+2.1725*+1.60089+1.17289+0.4320+0 .096 


H(z 


79:° 24437925 792? z 79 
6.351 £2 —5.99324+9.145 29 = 6.1152°43.3672— 1.029 


0.0902° +0 .5977°+1.3432*+1.792°+).34327+0.537240.09 


7.62026 —10.3382%+16.26724 = 14.10429 +9.8902* =4.524241.257 


0.04527 +0.3332°50.999:° +} 56524 +] 56572740 9392740 313240 045 





9.1642 = 16.225 2° 427.4552° — 28.8322 4+ 24.1092° — 14.5612°46.1342—1.522 


0.0222%40.1792"°+0 6272%4+1.254:5 4) 56824+1.2547°+0 6277740.1792+0 022 


1.8416 





0.033 
7 8! +0.923 89 + 2.17659 +1 429 84 + 1.358824 0.54982 +0.2148+0.031 








8°+0.9288°+1.9318*+1.2028° +0.9398° +0.3078+0.069 











3° +0.6968° +2.2428° 41 "2128° +1.5808*+0.596s540.35989 
+0.073s+0.013 









0.3687+0 06! 
6.376 2° — 10.0542 >415.73474=—14.6562° +10.6932°=—5.0858+1.502 


0.022 
8 39+1.1468  +2.6578°+2.1498°+2.1848°+1.14989+0.67482 11.0352 —24.106.2' + 44.4352 4 54.247 24 452.4624 = 398 .4342° +21.2647- 28.182 
+0.1539+0.024 
1-dB ripple (€ = 0.5088) 
N H(s H(z 
5 0.123 0.123:540.614:*+1.232541.23274+0.6142+0 123 
3°+0.9378*+1 689 8° +0.97482+0.5818+0.123 $.3032° —6.16924+8.9362°—6.6312°43.7252— 1.235 
6 0.06) 0 061 =°40.3687° 2441.228°°40.921:? 


27 915-4 75: 7S: - * 
7.6798: <15. tess 24222 — 29.121 24425.1267°=—15.81224 
+6.9192-1.816 






9.2542 —22.255:" +41.99326 —53.661 2° +53.5172*—40.6112° 
+23.24227=-9.2712+2.196 





H(z 
0.0827*+0 409244+0.817:°4+0 8177740 409 7+0.08} 
4.441 2° 6.1392 44+8.76329 = 7.07724 +4.1042—1.477 
0.0412°+0.2452°+0 6127%+0.8162°+0 6)2: 2+0.2457+0 041 
§.3472* =9 .6042°4+15.21024=15.0742° 411.2 297 24—5.677841.790 


0.02077 40 143:540 4282540 7147440 7142540 4282740 143740 020 
6.44525 14.2422 &+25.03425 —29.2032* +26.06 225° = 17 .0852447.7658—2.165 






0822" +0.2862°4+0.57}2°4+0 7147440 571 2° 40.2867? 
T7722 —20.3972 +39.5932°—52.7872° +54.37 12 * 42.6992 
+25.302? —10.46624+2.615 


0075 
8 8° 40.9208 '+2.42389+1.6668° +1 .8378*+0.8478° +0 4488 
+0.1078+0.017 
2 - dB ripple (¢ = 0.7648) 
Jee His 
5 0.08? 
9° +0.7078* +1.5002°+0.69484 +0.4598+0.082 
6 0.041 
8940.70 189 +1 746.84 +0.867 89 4+0.77284+0.2108+0.05) 
7 0.020 
81 40.6965 4+ 1.99485 +1 039 84+ 1.449° +0.38382+0.1672+0.020 
8 9.010 
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Step 2: Again, the design chart is based on a ripple edge frequency of 7/2, necessi- 
tating the normalization of the design specification frequencies. Looking 
at the Table of Digital Filter Frequency Transformations (Table 2.5), the 
following frequency transformations apply for a bandpass filter: 

Prototype filter from desired filter: 


j20, _ 2ak 76, | kad 
Peete Bales spur 


. 2.36) 
_ 2ak 276! k—1 326! ( é 
Pp” tae eect 


where: 


cos (6), /2 + 6,/2) 


~ cos (61, /2 — 6/2) 
Cad. 
k= tani") 2)icot (3 - 4) 
6, is the stopband frequency used to determine the filter order N. 
6. is the ripple edge frequency on the design chart, in this case 0, = 7/2. 


For this problem: 


Ca 
6 = 0.1335 
G, = 0.67 

; = 0.47 
§. = 0.57 


_ cos(0.37 +0.27) _ cos(0.57) 
~ cos(0.37 —0.27) cos (0.17) 
k = tan(7/4)cot(0.17) = 3.07 
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Since a = 0 


ite Oe te ee +0509 
1+ Et ei28 1 + (0.509)e2% 
522 Ue a 


~ 1+40.509e79-836 ~~ 1.394¢/9-282 


=> 96, = —2.86 rad = 2.86 rad (due to symmetry) 
Step 4: Determine the lowest order Chebyshev lowpass filter prototype from the 


design chart that gives: 
—3 db at 0, = 0.57 rad 
less than — 40 db at 0, = 2.86 rad 


Turning to Figure 2.11, the filter order is determined to be N = 3. Again, 
this agrees with the order obtained using the approach in the previous section. 
From Table 2.8 the lowpass prototype transfer function for N = 3 is: 


0.4912? + 1.4732? + 1.473z + 0.491 


Hupe(2) = 37723 — 1.oT72 4 2.047 — 0.750 


(2.38) 


Step 5: As with the Butterworth design, the Digital Filter Frequency Transfor- 
mation Table 2.5, is used to convert the lowpass prototype filter to a 


bandpass filter. 


Hgp(z) = Hzp,(z) Pe en 
oe See 


z= 
Zak ,4 b=1,2 
ripest easy 


_ 0.49123 + 1.47327 + 1.473z + 0.491 
Pole? — 1.2nree 4.9 2472 = 0.759 
_ 0.012z* — 0.0362* + 0.036z? + 0.012 


2.39 
z® + 2.13624 + 1.76822 + 0.540 ( ) 


Again, this is nearly identical to the transfer function obtained using the 


“traditional approach”. 
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Step 6: Verification of design using a computer generated frequency response 
plot. From Figure 2.13 it can be seen the design specifications of: 
C, = OAT = "ieto mrad 
@, = 0.67 = 1econad 
G, = 01337 — 041s rad 
are met. 
2. Elliptic Filters - A Direct Design Approach 

As was shown in Example 2.3, the design of elliptic filters is algebraically 
intensive, an extremely undesirable feature, especially for high-order filters. Elliptic 
filters characteristically exhibit a sharper cutoff, that is, a narrower transition width 
for a given filter order than their Butterworth or Chebyshev counterparts. For this 
reason they are very popular, making the pursuit of a more efficacious design 
approach a worthwhile endeavor. 

Just such an approach was found by taking the analog prototype filters of 
Table 2.4 and normalizing them to have a passband ripple edge (critical) frequency 
of one. They were then transformed to digital versions using the bilinear transfor- 
mation, resulting in digital prototype filters with a critical frequency of 7/2 (Table 
2.9). 

Thus, the amount of calculation involved in the design of a digital elliptic 
filter is cut in half because the requirement for finding an analog prototype is 
eliminated. The designer need only convert the filter design specifications to digital 
frequencies, find the digital lowpass prototype filter that meets these requirements 
(from Table 2.9), and then find the actual filter transfer function through the use 


of the digital frequency transformations of Table 2.5. 
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Two examples will be presented that illustrate this process; however, a 
more detailed explanation of how the transfer functions of Table 2.9 were arrived 
at is in order. 

For purposes of explanation, a second-order filter with 0.5 dB passband 
ripple, and a stopband gain of —20 dB will used. 


Turning to Table 2.4, the lowpass prototype transfer function is of the 


form: 

where, 
Ho = 0.100220 
Ao, = 5.33789 
Boi = 0.566660 
By, = 0.809390 

therefore, 


0.18? + 0.5338 


a See 
LPe(8) = 379 8004s + 0.5607 


(2.40) 


This prototype has a value for R of 2.76261. The parameter R is the ratio 
of the stopband frequency, w2,, to the passband cutoff frequency, w,,(i.e., R = 
w2,/ Wip). Thus, it describes the sharpness of the transition region. A large value 
of R is indicative of a broad transition region, while, conversely, a small value 
corresponds to a narrow transition region. As expected, the higher the filter order, 
the narrower the transition region, and the smaller the value for R. Figure 2.14 
(due to [5]), illustrates the magnitude squared frequency response of the normalized 
lowpass elliptic filters of Table 2.4. 

For this particular example, w;, and w2, are subject to the constraint 


that their ratio R be 2.76261. Furthermore, it should also be noted, the filters 
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transition 


Z 
lH(ja| passband band stopband 





Figure 2.14. Magnitude Squared Frequency Response 


of a Normalized Lowpass Elliptic Filter 


comprising this table have been normalized, so that the geometric mean of wi, 
and we, 18 one. 


(wiewae) = 1 (2.41) 


Combining these two constraints, the following relationships between w,, and wep 


apply: 


wip = 1/VR = 0.6016 rad/s (2.42) 


wep = VR = 1.6620 rad/s (2.43) 


The goal of the direct design procedure is to find a digital parallel to the 
analog prototype filters, and tabulate a table of digital lowpass prototypes. 

Step 1 of this process is to normalize the transfer functions of Table 2.4 
so that they all have a value for w;, of one. This is depicted Figure 2.15. 


For this example: 
0.13? + 0.5338 


H ee 
LPp(s) = 5 + 0.80943 + 0.5667 


(2.44) 


74 





Griginadl Normalized 


Figure 2.15. Normalization of Elliptic Filter 


where 


Wip = 1/VR = 1/V2.76261 = 0.6061 


To normalize H,p,(s) to wp = 1, let s =s/VR = 0.6061s 


0.1(0.6061s)? + 0.5338 
(0.60615)? + 0.8094(0.6061s) + 0.5667 
2 0.03675? + 0.5338 
~ 0.367452 + 0.49065 + 0.5667 


Pere ow, W), — 1, and wo, = A, since R = wo, /wi,- 


Hypp(s) = 
(2.45) 


The normalized analog prototype was then converted to a digital prototype 


through the use of bilinear transformation. 


Aypp(z) = Hrp,(s) 





+ 
" 0.0367 (51)° + 0.5338 (2,46) 


~ 0.367 (2=1)” + 0.4906 (251) + 0.5667 
_ 1.5752? — 2.752 + 1.575 
3.912? + 1.132 + 1.22 


The digital passband critical frequency is 6;, = 7/2, since an analog 





frequency of w1;, = 1 corresponds to digital frequency of 7/2 when using the 


bilinear transformation. 


Table 2.9 contains the compiled results of applying this process to some of 
the analog prototype transfer functions of Table 2.4. As with analog elliptic filters, 
special relationships exist between the passband critical frequency, 9,,, and the 
stopband frequency, 62,, of digital elliptic filters. 

Using the digital frequency transformations of Table 2.5, the following 


values for 0,, and 6), are found based on the analog prototype values for w;, and 


W2p- 
For 6;,, 
jo, _ JIt1 _ lei™/4 meee, 
= Se eee ee (2.47) 
= irs = n/2 
For 62,5, 
392 jk +1 Ter ak 1 j2tan7!R 
= OO EC ———— = 1 
—~jR+1  le-Jtan-1(-R) (2.48) 


Again, & is the transition region parameter described by the ratio of we, 
to w;, of the analog version of the filter. Table 2.9 presents a summary of the 


values for 6;,, 02,, and R for some of the digital lowpass prototype elliptic filters. 


Example 2.6 
Once again the design example of the previous section will be used to 


illustrate the new direct design procedure. The problem statement is: 


Design an elliptic lowpass filter that meets the following specifications: 


passband ripple of 0.5 dB 


passband ripple-edge frequency of 2 kHz 


e stopband gain of at most —20 dB for frequencies greater than 6 kHz 


sampling frequency is 20 kHz. 
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TABLE 2.9 
ELLIPTIC LOWPASS PROTOTYPE FILTERS 
Passband ripple = 0.5 dB ; Stopband gain = —20 dB 


(due to reference 5) 


Nip Gp A 

2 1.95708 2.447 2.7626 
3 1.5708 1.9157 1.9157 
4 1.5708 1.6145 1.0447 
0 1.0708 1.5862 1.0159 


Transfer Functions 


l= 


H(z) 
1.5752742.75z+1.575 
3.9 1s2 113s +122 


1.2762°+2.36727+2.376z+1.276 
4.626z2+0.008z2+3.030z—0.352 


1.4 179z4+2.366 z°43 .43622742.336z+1 -4179 
5.885 z4—1.09582z9+6 .5858z2—1.0954z+1.338 


1.7942°42.951224+4 85322944 .853227+2.9512z+1.794 
7.951625 —2.3142?4+12.719223—3.185622+4.9276z—0.902 


ao FF Ww Wd 
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Step 1: Find the critical digital design frequencies. 


Qrw, 2n(2x 10%) 





ra melee Seer err me = 
Qrw2  27(6 x 10°) 
, = = —————__ = 0.6 d 
22 me 30 x 103 ee 


Step 2: To determine the filter order of the prototype filter, find a and 42,, which 
are defined as follows from Table 2.5. 
_ sin (66/2 - 6,2) 
= Sin (0/246, 79) 
where 6, is the critical frequency of the prototype filter i.e., 0, = 7/2. 


Therefore, 
i sin (0.57 /2 — 0.27/2) 
~~ sin (0.54 /2 + 0.27/2) 


_ sin (0.7854 — 0.31416) _ 0.45399 (2.49) 
~ sin (0.7854 + 0.31416) 1.09956 
= 0.5095 


Find 62,, the stopband frequency for the lowpass prototype filter. 


wy) 
2p 


J92p = 
e 56 


1 — ael”2p 

34855 _ 0.5095 

~~ 1 = (0.5095 )e41-855 (2.50) 
= 1.255¢7°-789 1ei2-6785 

~ 1,.955e-I0.3965 


==> 05. = 20/50 rad 
Step 3: From the design chart (Table 2.9), determine the lowest order elliptic 
lowpass filter prototype. Table 2.9, indicates that a second-order fil- 
ter has a stopband frequency 62, of 2.447 rad, which meets the design 


specifications. Thus, the lowpass prototype transfer function is: 


1.5752 + 2.752 + 1.575 


H = 
LPp (2) = “sor + 1.13z + 1.22 


(2.51) 
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Step 4: Find the actual filter transfer function, using the prototype transfer func- 


tion and the digital filter frequency transformations of Table 2.5. 


Hy p(z) = Aipp(z) 





z—a 
l-—-az 


z= 








ae 
_ 1.5752? + 2.752 + 1.575 (2.52) 
Bee JL | see ay re: 250.5095 
™ T-(0.5095)z 
_ 0.58292? + 0.25412 + 0.5829 
~ 3.65122 — 3.80422 + 1.6593 
or 
0.1596z? + 0.0696z + 0.1596 
H,p(z) = —>_——____—_— (2753) 


z? — 1.0422 + 0.4545 
The transfer function of Equation (2.52) is identical to the transfer 
function obtained using the traditional design technique (see previous 
section). 

Step 5: Verify the design by obtaining a computer generated frequency response 
plot (see Figure 2.16). 

To further illustrate this direct design method, an additional example 
follows involving the design of a bandpass filter. 

Example 2.7 

A digital elliptic bandpass filter is to be designed to meet the following 
specifications: 
e 0.5 dB ripple in the frequency range 600 Hz < f < 900 Hz, 
e sampling frequency of 3 kHz, and 


@ maximum gain of -20 dB for 0 < f < 200 Hz. 


4, 
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3.14160 


N=2 


Step 1: Find the critical digital frequencies from the design specifications: 


ae ee 4 = 1.05 rad 
Oo TF, 3000 i al 
Qrf. 27x(900) 
a eae ig = 188 rad 
- 3000 : = 
Qn fo  2m(200) 
05. = oo = 3000 = (0.1337 = 0.418 rad 


6, = 6:0, = (1.25) (1.88) = 1.53 rad 


Step 2: Again, since the design chart is based on a cutoff frequency of 6. = 7/2, 
the design specifications need to be normalized to determine the lowpass 
prototype filter from Table 2.9. 

From the Table of Digital Frequency Transformations (Table 2.5), 
the following frequency transformations apply: 
Prototype filter from desired filter: 


j26 


, - of 
2 a 20k j9, k-1 
© ipa reo ee 


(2.53) 


elf2p — — =e 3 
oN 20k J 2 k-1 J 2 
reel Eom aa a 


where: 
_ cos(G,/2 + 4/2) 
~ cos (4, /2 — 6/2) 
k = tan (6~/2) cot (@,,/2 — 0, /2) 
2p 1s the stopband frequency used to determine the filter order, N. 


6, is the cutoff frequency of the prototype filter, 9, = 7/2. 
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For this problem: 


ys = 
6, = 0.1337 
2 — CU.en 

2 = 0.40 
0. = 0.57 


_ cos(0.37+0.27)  cos(0.5)a 
~ cos(0.3% — 0.27) ~— cos(0.17) 
a= 0 


k = tan(0.57/2) cot(0.6a/2 — 0.47/2) = 3.07 


Step 3: Find 62, to determine the lowest order prototype filter that may be used. 


Since a = 0 


Fe ris e!"2p + 0.509 
ee ie Sa 
1+ fete?» = 1 + (0.509)e??2r 
983640509 1.394e72:58 (ee) 


—j2.86 


= Ty (page = T gg4eIo7T = 
—— >, = 2-c0rat 
Step 4: From Table 2.9, determine the filter order of the prototype filter that 
corresponds to the above value of 92,. It can be seen from the table that 
a second-order filter fulfills the design specifications. 
Again, the lowpzss prototype transfer function is: 


1.57527 + 2.752 + 1.575 


H = 
LPe (2) = “397-34 1132 4 1.22 


(2.56) 


Step 5: As in previous design problems, the prototype filter transfer function is 
converted to the actual filter transfer function through the use of the 


digital frequency transformations of Table 2.5. 


pre) el jap (2) yoo oak hk = 
ae Sa ee SE 
aa Pes ee ee 


_ 1.5752? + 2.752 + 1.575 
36162 1132 = 1.92 eareniene 
-0.509z2-1 


_ 0.5833z4 — 0.25582" + 0.5833 
~ 3.650924 + 3.7996z2 + 1.6579 





2=—— 


(2.57) 





—_ 


Step 6: Obtain the frequency response plot for design verification (see Figure 
Del te). 

In summary then, this chapter has presented three commonly used recursive 
filter types: Butterworth, Chebyshev and elliptic. The type chosen by the designer 
is dependent on the requirements of his/her particular design problem. 

Butterworth filters exhibit a flat frequency response characteristic, but do not 
have as steep a transition region for the same order filter as do Chebyshev and 
elliptic filters. With Chebyshev and elliptic filters the steepness of the transition 
region is attained by accepting a degree of ripple in the passband. 

All three filter types may be designed using either the traditional approach, 
involving an analog prototype filter transfer function, or the direct design approach 


that uses digital prototype filter transfer functions. 


83 


4 


MAGNITUDI 
0.00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 14 


( 





~ —— 
ee 
a 
ae 
e 


| | ae | | 
| | : | 
0.00000 0.62832 1.25664 1.88496 2.51328 3.14160 


FREQUEN@ YH is 
N=4 


Figure 2.17. Frequency Response for 
Elliptic Bandpass Filter 


84 


II. NONRECURSIVE FILTER DESIGN 


A. INTRODUCTION 

The nonrecursive realization of digital filters is desirable because it can result 
in two very attractive features: linear phase and the absence of stability problems 
because all of the poles are at z = 0. 

Nonrecursive filters are normally characterized by a finite duration impulse 
response; consequently, historical methods of design involve the analytical deter- 
mination of the filter coefficients by expanding the desired frequency response in a 
Fourier series and then truncating the series to the desired filter length (order). A 
disadvantage inherent in this design method is the presence of an overshoot that 
occurs at discontinuities in the desired frequency response (known as Gibbs’ phe- 
nomenon). The use of window functions was devised as a remedy to this problem; 
examples are the Kaiser, Hamming and von Hann windows [14]. 

Another popular method of FIR filter design is frequency sampling, whereby 
the desired filter frequency response is sampled, and then the inverse discrete 
Fourier transform of these sample values is determined to find the filter impulse 
response. 

In this chapter both of these methods will be presented, and several design 
examples given, but first a review of the theoretical background of nonrecursive 


digital filter design by Fourier methods is in order. 
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B. BACKGROUND 
The design of nonrecursive digital filters hinges on the following relationships: 
A causal nonrecursive system can be described by a difference equation of the 


form [1] 
k=L 


y(n) = 5° bea(n — k). (3) 


k=0 


For an input z(n) = e?”® the steady - state system output is 
Yss(n) = €)"9 H(e?®) (322) 
where H(e?®) is the system’s frequency response. 


Expanding the right side of Equation (3.1) for the same input z(n) = e/”?, 


yields the following steady - state output 


Yss(7) = » beet ae = oy bp e387 eI *? 
k=0 k=0 


= bye?” + bye Fe F? + hy eIM e779 4 + bp cite IL9 a), 
=e”? [by + be 7? + bye 779 +... + bye FP9] . 
Equating Equations (3.2) and (3.3) gives 
eo)? (e3®) = 6)" Ib + bye 29 + bge 97? 4... 4 bye FE] (3.4) 
which implies 
H(e?*) = Ss bueno (si) 
But, by definition the frequency en an LTI system is 
H(e?°) = SS h(n)e~I"°, (3.6) 
n=00 


For a causal filter (h(n) = 0,n < 0), with a finite number of delays (n = L), 


however, the frequency response definition is 


n=L 
H(e!*) = S~ h(n)e7I"*. (35%) 


n=0 
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Comparing Equations (3.5) and (3.7) produces the following relationship that is 


the crux of nonrecursive filter design 
fe Ore (3.8) 


Equation (3.8) establishes a direct relationship between the unit impulse response 
of an FIR filter and the coefficients, b,, in the system difference equation. 

The goal of the design techniques presented in this chapter is to determine the 
filter coefficients (or weights), bo, 5,,..., 5, for a desired frequency response charac- 
teristic, Hp(e?°). These techniques include: Fourier Coefficient Design for lowpass, 


highpass, bandpass and bandstop filters; Windowing; and Frequency Sampling. 


C. FOURIER COEFFICIENT DESIGN 

As stated in the introduction we need to determine the filter coefficients, which, 
in the case of nonrecursive filters, are also the coefficients of the unit impulse 
response, h(n). Thus, a relationship between the desired frequency response and 
the impulse response needs to be established. 

Beginning with the desired frequency response 


Hp(e?’) = ys h(n)e?"?, (3.9) 


n= —oo 


It can be shown, [1], that h(n) of Equation (3.9) may be written, 
277+ ; ; 
h(n) = (1/2r) [ Hp(e}®)e"9d0. (3.10) 
9 


Equation (3.9) is recognized as the Fourier series expansion of the function Hp(e?°), 
where the A(n) are the Fourier coefficients (hence the source of the name for this 
design technique). 

Depending on the type of filter that is desired, this expression can be reduced 


to the following forms: 
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Lowpass Filters: 
hpp(n) = (K/xrn) sin(n@,) 


Ixy = desired magnitude in the passband 
(Say 
nia 0, oe 


6. = the desired cutoff frequency 


The number of coefficients 1s truncated to correspond to the desired filter order; 
as expected, the higher the order, the better the approximation to the desired 


frequency response. 





Figure 3.1. Ideal Lowpass Filter Frequency Response 


Highpass Filters: 
hu p(n) = hyp(n)(-1)" (3.12) 
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Figure 3.2. Ideal Highpass Filter Frequency Response 


Bandpass Filters: 








hap(n) = [2 cos(n)] hrp(n) (3.13) 
where, 
= se — ; 6, = upper passband frequency 
Gy = ite aa ; 6¢ = lower passband frequency 





Figure 3.3. Ideal Bandpass Filter Frequency Response 


Bandstop Filters: 
| has(0) = K — hgp(0) (3.14) 


hgas(n) = —hgp(n);n = £1, +2,... 


8&9 





Figure 3.4. Ideal Bandstop Filter Frequency Response 


As can be seen in the previous impulse response expressions for the various 
filter types, they are all based on a lowpass filter prototype, hyp(n). Thus, the 
design steps for the ideal lowpass prototype filter can be summarized as follows, 
bearing in mind that once the lowpass filter coefficients are determined they may 
be transformed to the coefficients for a highpass, bandpass, or bandstop filter, if 
desired. 

1. Lowpass Prototype Filter Design Procedure 
Step 1: Translate the design specifications to those of a lowpass prototype. De- 

termine the desired critical frequency, 9., and the passband magnitude, 
1g 


Step 2: Find the lowpass filter coefficients given by Equation (3.8), that is, 
hpp(n) = by =(N/an)sin(nO.); n= 0,+1,4+2,..., 42 (Jalen 


with L =(N —1)/2. 


Step 3: Shift hzp(n) to the right by L terms to make the filter causal. 
App(n) = ([A/a(n — L)]sin{(n —L)9,.]; n=0,1,2,... ,2L "(sale 
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Step 4: Transform the coefficients, hyp(n), to lowpass, highpass, bandpass or 
bandstop, as desired. Implement the design, and compare the frequency 
response obtained with the original specifications to ensure that these 


specifications are met. 


As an example, suppose a highpass filter with a passband of unity gain for 
frequencies greater than 10 kHz is desired. The system sampling frequency is 50 


kHz. 
Step 1: Cutoff frequency, 6, = wT = uel = ee = 0.4% 
Passband gain, K = 1.0 


Translating to lowpass prototype specifications: 


6, =7-O0,; where 6, is the cutoff frequency of the 
lowpass prototype and 6, is the cutoff 
frequency of the highpass filter 


=—6, =7-H =17-0.40 = 0.67 


Step 2: From Equation (3.11), the lowpass prototype filter coefficients are: 


Dees lsomn yy (ier) a — Oye lace 2. 10 


for a 21-coefficient filter, while the highpass filter coefficients are obtained 
from Equation (3.12). The resulting frequency response is illustrated in 


Figure 3.5. 


ol 


n hip hyp 

0 0.600 0.600 
a ea 0.303 -0.303 
aay -0.094 -0.094 
ceo -0.062 0.062 
sat! -0.076 0.076 
+5 0.0 0.0 
+6 -0.050 -0.050 
Sef 0.027 -0.027 
aie) 0.023 0.023 
+9 -0.034 0.034 
+10 0.0 0.0 


D. WINDOWS 
It has been shown that the expression for the desired frequency response, 
Hp(e?®), can be written in terms of a Fourier series. Equations (3.9) and (3.10) 


are rewritten below for convenience. 


Hp(e?*) = h(n)e~2”"? dé (3.9) 
where a is 
h(n) = (1/27) I Hp(e!®)ein dd. (3.10) 


In actual implementation, however, the infinite sum in Equation (3.9) must neces- 
sarily be truncated, since a filter cannot have an infinite number of delays. Also, 
at points of discontinuity in the ideal frequency response (where the magnitude 
abruptly changes from 1.0 to 0.0, or vice versa), the Fourier series approximation 
cannot match the desired frequency response exactly, even if an infinite number of 


terms were possible. There is always an overshoot of about nine percent near 
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Figure 3.5. Frequency Response for Fourier Coefficient Highpass Filter 
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discontinuities. This overshoot is the well known Gibbs’ phenomenon, and is alle- 
viated through the use of window functions. The filter weights, h(n), are modified 
using one of several available analytical expressions (Hamming, Von Hann, etc.) 
[14]. 


h(n) = h(n) - w(n) (3.159) 


h(n) = modified unit sample response 
h(n) = original unit sample response 
w(n) = window function 


The modified value of the desired frequency response is therefore 


n=L 


Hp(e?*) = \- h(n)e~i"9 
ane (3.18) 
= > h(n) + w(n)e?”? 
n=—L 


again where, 
h(n) = (1/27) iL Hp(e}®)ei"9 do 


which exhibits a gradual roll-off, rather than the steep slope characteristic of the 
ideal frequency response. Thus, the tradeoff involved when using windows in 
Fourier design is that a reduction in the overshoot caused by Gibbs’ phenomenon, 
is achieved at the expense of decreasing the sharpness of the frequency response 
cutoff. That is to say, the passband and stopband ripples are suppressed at the 
expense of a wider transition from passband to stopband; specifically, there is a 
less steep transition. Figures 3.6 and 3.7 illustrate these points; they depict an un- 
windowed and windowed lowpass filter design, respectively. For convenience, some 


of the more popular window functions are summarized [14]: 
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e Hamming Window 


w(n) = 0.54 + 0.46 cos(n7/L) 
| (3.19) 
(Da 1 oe bs eens ws 6, 


e Von Hann Window 


w(n) = 0.50 + 0.50 cos(nz/L) 
(3.20) 
Pe One eG 


e Kaiser Window 
(9 2)})1/2 
[9B (3.21) 
ss) sella se aed 


where Io(-) is the modified zeroth-order Bessel function, and @ is a constant that 
specifies a frequency response tradeoff between the peak height of the sidelobe 


ripples, and the width of the main lobe. 


E. DESIGN OF A DIFFERENTIATOR 
The Fourier series design procedure is easily applied to the design of a differ- 
entiator, which is often used in signal processing applications to track the rate of 


change of a signal. The desired analog frequency response for a differentiator is 
Hp(jw) = jw. (3.22) 
Its counterpart in the digital frequency domain is 
Hp(e?*) = j6/T; 6 = wT. (3.23) 
Figure 3.8 illustrates the ideal magnitude and phase characteristic. 
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- x2 


Figure 3.8. Ideal Differentiator Frequency Response Characteristic 
The Fourier series design approach requires the determination of the Fourier 
series that approximates this ideal frequency response characteristic. 
Substituting the desired frequency response (Eqn. 3.23) into Equation (3.10) 


yields 7 
h(n) = (1/27) [ (j0/T)e" dé 


r 


=(j/2rT) | 6e7"°dé. 


(3.24) 


Evaluating this integral, the following expression is obtained for the unit sample 


response of a differentiator. 


h(n) = (=1)P (nd), eee 
(3.25) 
AG) a 


9§ 


Figure 3.9 depicts the frequency response of a 20th-order differentiator whose co- 


efficients for T = 1 are as follows: 


n h(n) n h(n) 
ee 10 1 21 

1 1/9 1 
2  -1/8 tee 3 
SG 14 1/4 
4-1/6 15 -1/5 
5 1/5 16 1/6 
6-1/4 ily Pare 
7 SVE ig 1/8 
8-1/2 19 -1/9 
9 1 20 1/10 
10 0 


Also, included is a 20th-order differentiator with a Hamming window, which, as 
can be seen in Figure 3.10, is a very good approximation. 

In summary then, the Fourier series design technique is very effective, and 
can be used in many applications such as the design of a differentiator illustrated 
here. The use of window functions is necessary, however, to reduce the Gibbs’ 
phenomenon effect, which introduces a sacrifice in the designed filter’s transition 


region. 


F. FREQUENCY SAMPLING 
The window design technique introduced in the previous section has a distinct 
disadvantage in that the computation of Fourier series coefficients for the desired 
frequency response can be impossible when an analytical expression for the desired 
filter’s frequency response is unknown or extremely difficult to determine. 
Frequency sampling is a technique originally proposed by Gold and Jordan and 
later developed by Rabiner et al. as a solution to this problem [9]. This method 


has two distinct advantages: 
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e It is possible to design a filter that approximates any desired frequency domain 
specifications, without ever having to determine an analytical expression for 
the filter frequency response, H(e?°). 


e The filter coefficients, h(n), can be determined using the IDFT algorithm. 

Before proceeding with an explicit summary of the steps involved, and a de- 
tailed example, a discussion of the theoretical background will be presented. 

In many filter applications a sharp cutoff amplitude characteristic and linear 
phase are desired. To ensure that these objectives are met when designing a filter, 
requires consideration of the number N and location of the equispaced frequency 
samples. Discussion of these parameters relative to their impact on the filter’s 
frequency response involves consideration of four separate cases [9]. 


Case 1 N odd, frequency samples at 


6, = Izk/N,k =0,1,2,...,N —1. 


Figure 3.lla. Frequency Samples for N Odd 


Case 2 N even, frequency samples at 


6, =2nk/N,k =0,1,2,...,.N—1. 





Figure 3.11b. Frequency Samples for N Even 


Case 3 N odd, frequency samples at 


(k + 1/2) 


C—O V 


SO = 41, 





Figure 3.llc. Frequency Samples for N Odd 


Case 4 N even, frequency samples at 





Epo ee) 


6,. 
‘ N 


ee) eee. I ade 
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Figure 3.1ld. Frequency Samples for NV Even 

It should be noted that Cases 3 and 4 do not readily lend themselves to inverse 
DFT computation using the FFT algorithm, since the first frequency sample is not 
at 0 Hz. For this reason a detailed discussion of these two cases will be omitted, 
however, information regarding these cases can be found in [9]. A discussion of 
frequency sampling considerations for Cases 1 and 2 follows. 
Cases 1 & 2: Sampling Theorem 

Recall that for an FIR filter with impulse response, h(n) = {h(0), A(1),..., 
h(N —1)} , the transfer function is 


N-1 


IG] = 0 oe (3.26) 


n=0 
As discussed earlier, the impulse response, h(n), can also be represented in terms 
of its Discrete Fourier Transform (DFT), since it is of finite duration [9] 
N=-1 


h(n) = (1/N) }> Hyei? toi (32% 


k=0 


where Ay = H(z) 





zest 2ek/N 
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Substituting the expression for h(n) into Equation (3.26) gives [9] 


N-1 N-1 
H(z)= )> am metre z7 


n=0 k=0 
N-1 [N-1 
=(1/N) >_ bs Herel ao 
n=0 Lk=0 
ee (3.28) 
=(1/N) yy \- ge. a) ae 
k=0 n=0 
N-1 N-1 
=H) eee ae 26en 2-" 
k=0 n=0 


Applying the finite geometric sum property to the inner summation 





Nei 1 een)” 
. = z 
3 elarkn/N , Wes 
i _ ei 2xk/N 
n=0 : 
N-1 1 = | Sr | 
. nis Fe 
S | el2mkn/N , no _ 
1 = ei 2xkiN 
n=0 


i j2nk , —-N 
— | — citxk/N, 1 SS k= 0, i. Oe 


Substituting back into Equation (3.28) yields (9] 


| ‘Ele 
pn eee 


molec) _ — hey 


Let z = e/°, to determine the interpolated frequency response 





i8) eI Ne/2 (ei NO/2 _ e—Ind/2) iv Fy 
— (ej2mk/N\ (e—j8 
N aa (e/ Eien”) 
_ paTNO fe | ie. (ciN0/2 = e~J0N/2) 
_ — (ej2mk/N) (o—j8 
la (C7 Nea”) 
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oni N G2 ae (iN 0/2 _ ¢-iN0/2) . e-ink /N . 95/2 
e—jmk/N . eI 9/2 _ ei tk /N . e—J9/2 


ON 
k=0 
en3 8/2 eee Na) ae (ei NO/2 = gh”) seni mRIN 
~ N ej(9/2—ek/N) _ ej(0/2—nk/N) 
eWINO 60? SS nny , He sin (3) 
ye eb (3.29) 
sin (2 = a) 


N k=0 
Thus, Equation (3.29) expresses the interpolated frequency response in terms of 


the sample values, H;,, of the desired frequency response [9]. 
In Case 1, where the number of frequency samples, N, is odd, choosing the 





frequency samples, H;,, to be real and symmetric yields a real and symmetric 


impulse response 
(N-1)/2 
H 2 
=) ai cos (Fins) (3.30) 


h(n) = W _ N 
k=1 


The interpolated frequency response, H (e? ae is also purely real which, as stated 


in the beginning of this section, is highly desirable for most filter applications 
(N-1)/2 
(3.31) 


H (e") = > 2h (n)cos(n@). 
no et 


Case 2, however, presents a problem. Looking at Equation (3.29), it can 
~itk/N inside the summation introduces an 


be seen that for N even, the term e 
imaginary component into the interpolated frequency response, H (e? Ne 


For example, suppose the following lowpass filter frequency response is sampled 


from —x < 0 < =z with an even number of frequency samples, N = 4. From 


Equation (3.27), the following noncausal impulse response is obtained 


nay = (174 Hypeune 
)=(1/4) S> A aan 


k=-—2 
= (1/4) fo BA( len! ompeeale ie! aan 
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Figure 3.12. Desired Lowpass Filter Frequency Response 


therefore, . 
h(—2) = (1/4)[(0 + &* +1 +e777] = —0.25 


h(—1) = (1/4)[0 + ef"? +14 €797/?] = 0.25 
h(O) = (1/4)[0 + 1+1+4 1] =0.75 
h(1) = (1/4)[0 + e72*/?7 4.14 7/7] = 0.25 


and the frequency response is 


me y= >) A(n)je**? 


n=—2 
= (—0.25)e??? + (0.25)e7® + (0.75) + (0.25)e77? 
= (—0.25)[cos(2¢) + j sin(26)] + 0.25[cos(6) + j sin()] + (0.75) ‘238) 
+ (0.25)[cos(6@) — 7 sin(@)] 


= (0.75) + (0.5) cos 6 + (—0.25)[cos(2) + 7 sin(24)). 
The imaginary component introduced is 
—0.257 sin(26). 


According to Reference 9, the amplitude of this component should be, 


A=(1/N) }> Hi(-1)' 


k= —2 
= (1/4)(0 —1+1-—1] = -0.25 
which is indeed the case. 


If an odd number of samples, N = 5, were taken of the same frequency re- 


sponse, the impulse response would be 


ha) = (yee eee 


k=-—2 





ern ay fs 3 3754 
= (1/5)[0 + (1)e7? +1 + (1)e 75" + 0] Geo 
= (1/5)[{1 + 2cos(27n/5)| 

therefore, 
0) Ss 

h(+1) = (1/5)[1 + 2 cos 27/5] = 0.324 

h(+2) = (1/5)(1 + 2 cos 4x/5] = 0.124. 
The frequency response is thus 

. - ° 
He?) = Giese 
n=—2 
= (0.324)e7?? + (0.124)e7° + 0.6 + (0.124)e72 + (0.324)e77?° 
(3.35) 


= (0.324)[cos(26) + j sin(20)] + (0.124)[cos 6 + 7 sin 6]+ 
0.6 + (0.124)[cos 8 — 7 sin 4] + (0.324)[cos(20) — 7 sin(24)] 
= 0.6 + (0.648) cos 0 + (0.248) cos(26) 

and there is no imaginary componc.it. 


To remedy the imaginary component problem for the even case, the following 


substitution is made for the frequency sample values Hy [9] 
H, = Gyel™*/%, (3.36) 


The set of G,’s is selected so that G(0) = H(0), G(N/2) = 0, 
and G(k) = —G(N —k), k =1, 2,..., (N/2)—-1. 


The resultant impulse response is both real and symmetric, with h(n) =h 
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(N —1—n);n=0, 1, 2,..., (N/2) —1. The values of h(n) are given by, 


h(n) = a FE oo | (5 Jk + (FF) ma (3.37) 


k=1 
and the interpolated frequency response H(e?®) is, 


(N~1)/2 


H(e?#)= 5° — 2h(n)cos(n8). (3.38) 


n= — (N=) 


The steps involved, when designing a filter using the frequency sampling 
method, are as follows: 
Step 1: Given the continuous frequency response specifications for a desired filter, 
N samples are taken at equispaced frequencies over one period of the 
desired response 


H(k) = H(e?*); 6 = 2nk/N. (3.39) 


If the number of samples, N, is even, they need to be transformed using 


the following relationship, prior to proceeding with Step 2 


H(k) = G(k)e??™PF/N (3.40) 
such that 
G(0) = 4H(0) 
G(N/2) =0 
G(k) =-—G(N -—k);k =1,2,...,(N/2) —1. 


Step 2: Determine the Fourier coefficients of the desired filter by finding the 
IDFT of these sample values. In other words, determine the impulse 


response. 


Ga sCN pe riche 22257 (3.41) 
N 
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Step 3: 


Step 4: 


Step 5: 


Again, it may be necessary to use an appropriate window function, 
h(n) = h(n) - w(n). (3.42) 


The Fourier coefficients thus determined correspond to the weighted filter 
impulse response, h(n), and the filter is realized nonrecursively by the 


difference equation 
N-1 
Gee ye »y b,x(n — k) (3.43) 
k=0 


where 


by = h(n). 


Verification of the filter design is accomplished by determining the inter- 
polated frequency response, H(e?°), resulting from the use of the above 
filter coefficients. This frequency response is compared to the original 
desired frequency response, Hp(e?®), to see if it is a reasonable approxi- 


mation. 


This procedure is illustrated by the design of a bandpass filter. First, an odd 


number of frequency samples will be used and then an even number. 


Example 3.1 Frequency Sampling for N Odd 


A bandpass filter is desired with the ideal frequency response characteristic 


illustrated below. The number of frequency samples to be used is N = 255. 


Sampling gives 


Oy = 2rk/N 


where N = 255 
(3.44) 
b= Ole ee oe 


or AO = 2n/N = 0.0246. 


110 





Figure 3.13. Desired Bandpass Filter Frequency Response 


Thus, 
Ok 62 OO ena k= 149. 2... 191 
H(s) = | 00 k= 054. 60 ane ee VON els 
and A= 192.,.<, 204 . 


Taking the IDFT of the above frequency sample values and truncating h(n) to 51 
terms with a rectangular window, produces the interpolated frequency response of 
Figure 3.l4a, while Figure 3.14b illustrates the response obtained using a Hamming 
window. 
Example 3.2 Frequency Sampling for NV Even 

Repeat the previous example, with N = 256 frequency samples. 


Sampling gives 


0, = 2nk/N 
where 

N = 206 

kK = 50, lee 2 00 
or 
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Figure 3.14a. Unwindowed Bandpass Filter Frequency Response 
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Figure 3.14b. Windowed Bandpass Filter Frequency Response 
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Thus, 
1.0; k= 64,2... 106 sand (= 1507 
H(s) = {00 k= 0203 and k = Wi....,14e 
andi k = 193... 200". 


But recall for N even, the sample values H(k), must be transformed to elimi- 


nate the unwanted imaginary component. Applying Equation (3.40) gives 
H(k) as Gikje2 es 
such that, 
C(O) =O) = 0 
G(256/2) = G(123)} =0 


G(k) = —G(256 — k);k =1,2,...,127. 


For example 
G(65) = H(65) = 1.0 = —G(191) => G91) = —10 


and the transformed H(65) and H(191) are 

H(65) == 1.0¢e/27(65)/256 

H(191) = —1.0e)7*02256 
Using the transformed H(k) values, the IDFT is determined to find the impulse 
response, which is again truncated to 51 terms. Figures 3.15a and 3.15b show the 
unwindowed and windowed frequency responses, respectively. It should be noted 


that when using a Hamming window for an even number of terms (in this case 52) 


the equation for the window function (Equation 3.19) should be modified as follows 


N-1 


w(n) = 0.54 + 0.46 cos ome] (3.45) 
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Figure 3.15a. Unwindowed Bandpass Filter Frequency Response 
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Figure 3.15b. Windowed Bandpass Filter Frequency Response 
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G. TRANSITION POINTS 

As discussed in the previous section, frequency sampling involves taking N 
equispaced samples of a desired filter frequency response, Hp(e??). The IDFT 
algorithm is used to approximate the unit sample response, h(n), which, for the 
case of nonrecursive filters, is equivalent to the Fourier series coefficients. The 
coefficients, in turn, are used to design a filter that approximates the original 
desired continuous frequency response. 

A filter designed in this manner has a frequency response that is equal to 
the desired frequency response at the frequencies of the sampled values, however, 
the response can deviate significantly from the desired response at frequencies in 
between these sample values. 

A method that can be used to smooth the frequency response involves the 
use of transition samples between the passband and stopband of the desired fre- 
quency response [9]. The values of these transition samples are selected based on 
minimization of the ripple in the passband and/or stopband, or minimization of 
the maximum sidelobe of the frequency response. In other words, a minimization 
algorithm is used to find the optimum values for the transition samples based on 
the selected minimization criterion. 

As before, an example will be given to illustrate this technique, but first a 
summary of the steps involved: 

Step 1: Sample the desired continuous frequency response at N equispaced fre- 
quencies where JN is the number of Fourier coefficients that will be used 
to approximate the desired filter response. The spacing between the 
frequencies is A6 = 27/N. 

Step 2: Determine the impulse response h(n) (Fourier coefficients) by finding the 


Inverse Discrete Fourier Transform, IDFT, of the sample values. 
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Step 3: Using the coefficient values determined in the previous step, find the con- 
tinuous frequency response and compare with the desired filter response. 

Step 4: Use the minimization algorithm to adjust the transition coefficient values, 
in order to obtain as close a match as possible between the desired filter 
response and that obtained through the design procedure. 

The minimization algorithm as used in [9], is based on minimizing the maxi- 
mum sidelobe of the frequency response. Tables were generated for lowpass filters, 
bandpass filters and wide-band differentiators of varying bandwidths; 1, 2, or 3 
transition points; and odd and even values of the number of frequency samples, NV. 

Table 3.1 is a reprint of subsets of these tables (due to Reference 9) for lowpass 
filter design, using 1, 2, or 3 transition points. The number of frequency samples is 
N = 15, the column labeled minimax refers to the maximum sidelobe, and Type-1 
data means that the first frequency sample is taken at 0 = 0. 

Design Example 

In this example the goal is to design a lowpass filter with three frequency 
samples in the passband symmetric about the origin (BW = 3), and a total of 15 
frequency samples (N = 15). The effects of using 0, 1, 2, and 3 transition points 
will be investigated. 

The values of the frequency samples, transition point values, and impulse 
responses, h(n), are summarized below for all four cases; the transition point values 


were obtained from Table 3.1. 
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BW 


me © WH Fe 


TABLE 3.1 
LOWPASS FILTER DESIGN 
(TYPE-1 DATA, N ODD) 


ONE TRANSITION COEFFICIENT 


(due to reference [9]) 


BW 


ao Fr |W NO Fe 


Minimax 


N=15 


-42.30932283 
-41.26299286 
-41.25333786 
-41.94907713 
-44.37124538 


6 -56.01416588 
TWO TRANSITION COEFFICIENTS 


Ti 


0.43378296 
0.41793823 
0.41047363 
0.40405884 
0.39268189 
0.35766525 


BW 


mm GW ho Fe 


3) 


Minimax 


Ty 


[il 1s 


-70.60540585 
-69.26168156 
-69.91973495 
-75.51172256 
-103.46078300 


0.09500122 
0.10319824 
0.10083618 
0.08407593 
0.05180206 


TL? 


0.58995418 
0.59347118 
0.58594327 
0.55715312 
0.49917424 


THREE TRANSITION COEFFICIENTS 


Minimax 


-94.61166191 
-104.99813080 
-114.90719318 
-157.29257584 


Ty 


NGS 15 
0.01455078 
0.01000977 
0.00873413 
0.00378799 


Jhis) 


Tp 


0.18457882 
0.17360713 
0.16397310 
0.12393963 


T3 


0.66897613 
0.65951526 
0.64711264 
0.60181154 


Case 1 - 0 Transition Points 
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Figure 3.16a. Case 1 - Frequency Samples 


ku 20 +1 +2 +3 +4 +5 +6 i7 
H(k) 3 1 1 1 0 0 0 0 0 
h(n) — 0.333 0.278 “0.142 Ol —0077 e=Or0o7 10 OmOrn5e 


Case 2 - 1 Transition Point 


N= 15 
T,= 0.410474 
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Figure 3.16b. Case 2 - Frequency Samples 


k — 0 acl 2 +3 +4 +5 6 aay 


H(k) > 1 1 1 0.410 0 0 0 0 
h(n) — 0.388 0.295 0.098 -0.044 -0.061 -0.012 0.017 0.014 


Case 3 - 2 Transition Points 


13 
0.100836 
>= 0.585943 


N 
i 
a 





Figure 3.16c. Case 3 - Frequency Samples 


ie Sen0 +1 42 43 +4 +5 +6 +7 
H(k) > 1 1 1 0.586 0.101 0 0 0 
h(n) — 0.425 0.300 0.066 -0.059 -0.041 0.005 0.013 0.004 


Case 4 - 3 Transition Points 


= 15 

.008734 
nO S 973 
= 0.647113 


ll 
© 


N 
T,= 0 
aR 
ab 





-7 -6 -$ -4 -3 -2 -1 Meco wetmoo? 8 (k 
k — 0 aenh a2 ao +4 EO ==0 Sef) 
H(k) — 1 i 1 0.647 0.164 0.009 0O 0 


h(n) — 0.443 0.301 0.050 -—0.062 -0.032 0.008 0.010 0.002 


Figure 3.16d. Case 4 - Frequency Samples 


Figure 3.17 depicts the results of frequency sampling design using transi- 
tion points. As expected, the amplitude of the passband and stopband ripple is 
reduced with an increasing number of transition points; however, the tradeoff is 
that the sharpness of the cutoff is reduced. This example illustrates the usefulness 
of minimization algorithms in smoothing interpolated filter frequency responses 
obtained with the frequency sampling technique. As stated earlier, Reference 9 
contains tables for lowpass and bandpass filter design, as well as wideband differ- 
entiators. While quite extensive, a designer may want parameters that have not 
been tabulated. In this case, approximate values of the transition coefficients can 


be derived from the tabulated values given, using linear interpolation. 


H. DESIGN OF AN INTEGRATOR 


In Section A analytical methods were used to design a differentiator. The use of 
these methods to design an integrator, however, does not yield an easily obtainable 
solution. As will be shown shortly, frequency sampling solves this dilemma by 


providing a straightforward design technique that produces extremely good results. 


Bandpass Integrator Design Example 
A bandpass integrator with the following frequency response characteristic 
is desired. Both of the aforementioned design methods will be applied to the design 
problem to illustrate the advantage of the frequency sampling technique. 
a. Method 1 - Analytical 
In the s-domain an integrator is described by the transfer function 
H(s) = 1/s. The frequency response is thus: H(jw) = 1/jw = (jw)7'. Converting 


the frequency response to digital frequencies 


H(jw) =1/jw => H(e?*) = —jT/0;T = period. (3.46) 
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Figure 3.18. Ideal Integrator Frequency Response Characteristic 


To determine the sa as h(n), 


* * 62 . . 
h(n) =(1/2e) [Hp (3%) ei d0 + uh Hp (€2%) edo 
8 


62 i 


* @2 . 
kn) = ane [ " (=jT/0) 8" d0 + Bs 1/64 dé 
—82 A, 
Ve az az’ ax 
eas |B: MN as 


(jn)?6? (nye + —O, 





eS -jT/2z {In om 7 7 Ee Dino! 3-3! 





—65 
: : 8 
—jT CE eG ; 
(a) ec alee 


A considerable number of steps later yields 


h(n) =T/x nC —6;)+ is ET (he OP) 4+ ...4 5 (05 — a) (3.47) 


where n is the number of coefficients desired for the implementation of the bandpass 


integrator. 


The analytical expression for the impulse response of the bandpass 


integrator, Equation (3.40), is cumbersome, especially when the integrator design 


is of a high order. Therefore, an alternative method is desirable. 
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b. Method 2 - Frequency Sampling 


Recall that the expression for the desired frequency response is 
Hp (e?°) = -jT/8. (3.48) 


Letting T = 1, one hundred and one frequency samples will be taken of the de- 
sired bandpass integrator with a bandwidth of 0.17 to 1.67. In other words, 101 
frequency samples, H(k), of the magnitude and phase characteristic (Figure 3.19) 


will be taken using the following increment of the digital frequency 0 


NG = 2a = 0,002k for V = 101 


3.49 
H(k) = H (e*) Eye 506s 60. ee) 
6=0.062k 


The 101-coefficient impulse response is obtained by determining the 
IDFT of these frequency sample values, and a 21-coefficient rectangular window is 


applied yielding the following bandpass integrator coefficients. 


n  A(n) n h(n) 

0 0.0974 11 0.4813 
1 0.0793 12a 02239) 
2 0.0946 13 0.2248 
3 0.0439 14 0.0818 
4 0.0362 15 0.0562 
5 —0.0562 16 —0.0362 
6  -0.0818 04 39 
(0 2248 18 -—0.0946 
& Sih peel Lou (0793 
o) e=04813 20 —-0.0974 
10 0.000 


Finally, the interpolated frequency response is obtained. As can be 
seen in Figure 3.20, the results are good. (Solid dots indicate the ideal frequency 
response.) Also included in Figure 3.21 are the results obtained when a Hamming 


window is used, while Figures 3.22 and 3.23 are the results of using 41 coefficients. 


PST 





Figure 3.19. Ideal Bandpass Integrator Frequency Response Characteristic 


This example further demonstrates the versatility and usefulness of 
the frequency sampling technique. It allows for the design of filters that may be 


unattainable using traditional analytic methods, and is thus an indispensable tool. 
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Figure 3.20. Unwindowed 21-Coefficient Bandpass Integrators 


(Solid dots indicate ideal frequency response) 


130 


3.14160 


MAGNITUDE 


Zo 


4 \ 
3 ae 
2/ 


0.00000 0.62832 1.25664 1.88496 2.51328 3.14160 
PRONG Ya bis 


Figure 3.21. Windowed 21-Coefficient Bandpass Integrator 


(Solid dots indicate ideal frequency response) 
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Figure 3.22. Unwindowed 41-Coefficient Bandpass Integrator 


(Solid dots indicate ideal frequency response) 
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Figure 3.23. Windowed 41-Coefficient Bandpass Integrator 


(Solid dots indicate ideal frequency response) 
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IV. COMPUTER-AIDED DESIGN 


The basis of all computer-aided design (CAD) techniques is optimization. This 
is to say, a desired filter frequency response is approximated by a particular filter 
whose coefficients are to be determined. The accuracy of the approximation is 
evaluated according to some criterion, usually an error function, that indicates 
how large a disparity exists between the desired filter frequency response and the 
approximated filter frequency response. Variable parameters of the approximating 


function are then “adjusted” to optimize the filter design in terms of this criterion. 


A. REMEZ EXCHANGE ALGORITHM 

An extensively used computer-aided design technique for linear phase FIR 
filters is the Remez exchange algorithm. The mathematical basis for this algorithm 
is the weighted Chebyshev approximation. 

A summary of the approximating and error functions for this algorithm follows. 
It has been shown in [1] and [13], that the frequency responses for the four cases 
of linear phase filters - i.e., even or odd symmetry with an even or odd number of 


terms - can be written in the form: 
ie (€?°) = oe IHN-1)/2 04/2) L FF (€?°) (4.1) 


where H (e?°) is a real-valued function used to approximate the desired filter’s 
magnitude specifications, and the remaining terms approximate the desired phase. 


Table 4.1 gives the values L, and the form of H (e7°) for all four cases. 
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TABLE 4.1 
FREQUENCY RESPONSES FOR LINEAR PHASE FILTERS 


(due to reference 13) 





L ahife®) 
Case 1 - N odd 
(N~1)/2 
Symmetrical 0 S> a(n) cos(né) 
n=0 
Impulse Response 
Case 2 - N even 
N/2 
Symmetrical 0 5° &n)cos [6(n — 1/2)] 
n=] 
Impulse Response 
Case 3 - N odd 
(N-1)/2 
Anti-symmetric 1 S> c(n)sin(n6) 
n=1 


Impulse Response 


Case 4 - N even 
N/2 

Anti-symmetric 1 > d(n)sin[6(n — 1/2)| 
n=1 

Impulse Response 


Using trigonometric identities, the expressions for H (e? @) in Table 4.1 can be 


rewritten in the following form: 


H (ec?) = Q (e*) P(e?) (4.2) 
where Q (€7°) is a fixed function of frequency, 8, and P (e? ss consists of a sum of 


weighted cosine terms (see Table 4.2). 
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TABLE 4.2 


FREQUENCY RESPONSES FOR LINEAR PHASE FILTERS 


where 


(due to reference 13) 








Q (e’") P (e’") 
(N-1)/2 
Case 1 1 >> = @(n) cos(né) 
(N/2)=3 


Case 2 cos(@/2) 22, b(n.) cos(n8) 


(N-1)/2 

Case 3 sin(@) Ze, €(n) cos(né) 
(N/2)=1 

Case 4 sin(6/2) = d(n) cos(n@) 


ee [a(S 

a7 { evan heal for n = 1,2,..., 1 
b(1) = 6(0) + 1/20(1) 

b(n) so) = be 2)- fot i =)293 eel 
b(N/2) = 1/26 (¥% —1) 
c(1) = &(0) — 1/2é(2) 

. A for n = 2,3,..., % 4 
<(MaL) = sy 26( Sa 
d(1) = d(0) — 1/2d(1) 

d(n) a) = at fon ae eel 
d(F) =1/2d(F -1) 
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Since Q (e/®) is a function of frequency only, it can be seen that the approx- 
imating function, H (e?®), can only be optimized in terms of P (e/®); that is, the 
dependency of H (ei *) on the filter coefficients is contained in P (e/ °). Thus, the 
coefficients of P (e/ ¢ ) can be varied to achieve an optimum filter design, and, the 
true approximating function can be generalized for all cases as, 


L 


Ip ( eh) = = a(n) cos(n@) (4.3) 


n=0 

where L = (N —1)/2 or (N/2) — 1 depending on phicht eases being considered, 

and the a(n) are the weights from which the filter coefficients can be determined. 
As stated in the introduction to this chapter, a measure of how well the de- 

signed filter frequency response approximates the desired filter frequency response 

is required. The weighted Chebyshev approximation uses an error function defined 

as follows: 


E(8) = W(8) [Hp (e*) — # (e*)| (4.4) 


where 
Hp (e? "1 = the desired frequency response 


H (e?°) = the designed frequency response 
W(6) = weighting factor 


E(@) = error 
In order to see the relationship between the filter coefficients and the error, 


Equation (4.4) can be rewritten in terms of the functions Q (e/®) and P (e?%). 


E(8) = W(9) [Hp (7°) — Q (e?*) P (e?*)] 





= iy |Ho le") ay is (4.5) 
= W(8)Q (e*) | Deisy ~ P(e’) 
Letting 
W(6) = W(8)Q (e?°) and Hp (e?°) = Hp (e?®) /Q (e?°) (4.6) 
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yields 
E(8) = W(0) [Ap (c%*) — P (e”*)| (4.7) 


Thus, the Chebyshev approximation that is performed using the Remez exchange 
algorithm can be stated as follows: 

Find the set of filter coefficients (determined from the values of a(n)) that 
minimizes the maximum absolute value of the error, E(6), over the frequency range 


of interest. 


meen == mi | 0 =| (4 8) 


(coefficients ) 


At this point, a discussion of the weighting function, W(@), is in order. The 
purpose of this weighting function is to ensure a small tolerance for error in critical 
frequency ranges. 

If W(@) is large, this means a large deviation between the desired frequency 
response, Hp (e oN and the designed frequency response, H (e ay cannot be toler- 
ated. Looking at Equation (4.4) we see that if W(@) is large, the difference between 
the desired and designed frequency responses, lH D (e/ 2) =i (e *)| has to be small 
to keep the weighted error small. 

Conversely, if W(8) is small, the difference lH p (e??) - H (e7 ®)| can be larger 
and still meet the error criterion. Small values for W(4) would be used in frequency 
bands where close approximation to the desired frequency response is not critical. 

A copy of a program employing the Remez exchange algorithm [16] tas been 
included, and its use is best illustrated with an example. Before proceeding with 
the example, the salient features of the FIR Linear Phase Filter Design program 


can be summarized as follows. 
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1. Main Program Functions 
e Manage the Input - 


NFILT = filter length in terms of samples, 3 < NFILT < NFMAX 
NFMAX = 128, but can made greater or smaller by the user 
JTYPE = filter type 
1 = Multi-Passband/Stopband 
2 = Differentiator 
3 = Hilbert Transformer 
EDGE = number of frequency bands, specified by 8jower and OBupper 
( maximum allowed is 10 ) 
FX = desired frequency response magnitude in each band, |Hp (e?°) | 
WTX = positive weight function, W(@), in each band 
LGRID = Grid Density (the grid density, default value is 16) 


e Set-up appropriate approximation problem, based on the desired frequency 
response, Hp (e?°), and the weights, W(@). 


e Yield Output - 


o The coefficients of the best impulse response obtained from the best cosine 
approximation. 


o The optimal error (sin "3" (6) 
o The extremal frequencies where E(@) =o |E(6)| 


2. Most Important Subroutines 
EFF - defines the desired frequency response, Hp (e a 
WATE - defines the weight function, W(@) 
REMEZ - calculates the best Chebyshev approximation of the desired frequency 


response 
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3. Important Note 


It should be noted that the program output is the impulse response, h(n). 
The user must supply his/her own program to calculate and plot the filter frequency 
response, which is necessary to determine the suitability of the filter design. 
4. Input Format 
NFILT JTYPE # Bands LGRID 
EDGE (Band Edges) 
FX (Magnitude in Bands) 
WTX (Weights in Bands) 
5. Bandpass Filter Design Example 
We wish to design a bandpass filter that meets the following specifications: 


Passband with a gain of one for the frequency interval 10 kHz to 15 kHz, 


with a system sampling frequency of 50 kHz. 
e The digital frequency band is: 


_ 2 (10') 4 
> Saoay eee 


_ 2n(1.5) (10*) _ 
us 5104) = 0.57 


e Normalizing so that 6 = m corresponds to a normalized frequency of 0.5. 


0.47 Oe 

—$ = y= (0 
wT 0 “Sr ae 

0.67 0. 

— = — = 6, = 0.3 
wT 0.5 


140 


e Designate Band Edges, Magnitudes, Weights: (the weights were arbitrar- 
ily made equal for this example) 


Band 1 

Edges - 0.0 to 0.19 
Magnitude - 0.0 
Weight - 10.0 


Band 2 

Edges - 0.2 to 0.3 
Magnitude - 1.0 
Weight - 10.0 


Band 3 

Edges - 0.31 to 0.5 
Magnitude - 0.0 
Weight - 10.0 


e Select Filter Length ( number of coefficients ) 
ve 
6. Sample Input File 
00021 00001 00003 0 
mo 0:19 0.2 0.3 Band Edges 


0.31 0.5 
0.0 10 860.0 Magnitudes 
10° 10-0 =—:10.0 Weights 


141 


es 


The output file below was generated using the above input file. 


Sample Output File 


N dte 


the output consists of the filter impulse response, the input values used and the 


deviation between the desired and approximated frequency responses. 


SVT eo <P NG = ee | 
Reet TYDVIJT ST IEQDONST kee 
35 al =" O21599373e.09) =m 
Hf >} = Jeo 7232 15-0 aa 
4( 3) = 91199799 te oa — ae 
| i = =e Oe = 
H( 5) = =I, 18S5h53NTF4ND = Ht ( 
Si & = =. 273 4212s Oa a= 
nt | = VLBA ZBARUTS4+IA = uf 
HC 2) = =001 22) 3048 2B = 
H( 9) = =9.23 07990 H 
(19) = —-). F6ETFAINM=-RNG = 
A(11) = 9. 2831 ny 
Bo NDe24 Bae 
TOUT? SAN Cea ae e200 200) V2 9 
TPPeO RAND FOGF aaa 4. NSS 
UE Sy eee be Ore ov okelare vel Tee 
W2CG Hees Oa 1 Oe) 
DEY Lae nN O23 305s ae Oe 3k 
DEVIN TION EN as = 9 625 ho a7. Drea 
ENS cen ean ECU ENG 5255s op ape? anc 
0. USO e O22 Seta J. 156274Pr9 
V9 5 Tass A299 09 OS ye an 
Je FHC 2920 9.5999909 
8. Results 


SET to> THOR oo  f Oe aac 
LINE\? DHS Se) et lee 
PRM" SYCHVNG" NeESOe Le 


PNTEDOONS Soe ie ee 


(ero 


25 TGN 
*«&* 

21) 

20) 

1a) 

18) 

7 

14 

15 

14 

ae 

12 

we 
Bae BA 0 ae 
O00 4.31950 05 
aAngn Aan 
9990 0. 7009009 
AnAN 1.990999 
eer alg 2 3u4u6794 
eS ~9, 25 16990 
ne cyrys 
ae 1 MORO 0.29900) 
als BAM Vac) os 1. 39522 


Figures 4.1 to 4.3, illustrate the frequency responses for filter lengths N = 


21,41 and 61, respectively. 


As expected, the frequency response improves with an increasing number 


of coefficients. The cutoff is sharper and the magnitude more closely approximates 


0 dB in the passband. As with other design methods, windows may also be used 
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to improve the frequency response characteristic. Figures 4.4 to 4.6, show the effect 


of changing the band edges as follows, for N = 21 coefficients. 


Figure 4.4 - Band 1 
Band 2 
Band 3 


Figure 4.5 - Band 1 
Band 2 
Band 3 


Figure 4.6 - Band 1 
Band 2 
Band 3 


0.0 to 0.19 
0.2 to 0.3 
0.31 to 0.5 


0.0 to 0.15 
82 to. 0.3 
0.35 to 0.5 


0.0 to 0.1 
(2 to0:3 
0.4 to 0.5 


Figure 4.4 represents the original cutoff frequencies, while Figures 4.5 and 


4.6 show the results of not selecting the stopband frequency cutoff close enough to 


the passband frequency starting point, i.e., the passband is broadened beyond the 


desired 0.2 to 0.3 specifications. 


Thus, when selecting the stopband cutoff frequencies one should choose 


values as close to the passband frequencies as possible. It is of importance to note 


that the program does not allow for the upper stopband frequency equaling the 


lower passband frequency. In other words the following parameters would not be 


acceptable: 


Band 1 0.0 to 0.2 


Band 2 0.2 to 0.3 


Band 3 0.3 to 0.5 
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(DB) 
-20.0 -10.0 


MAGNITUDE: 
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10.0 
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9.00000 0.62832 1.25664 1.88496 2.51328 3.14160 
FREQUENCY, THES 


Figure 4.1. 21-Coefficient Remez Bandpass Filter Design 
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Figure 4.2. 41-Coefficient Remez Bandpass Filter Design 


145 


1.88496 2.51328 3.14160 
PREQUENG ] THETA 


1.25664 


0.62832 





0.00000 


0;OT O°'O O'OT— 0'08- 0'0F-- 0'OF— 0'0SG— 0°0980'02- 
(10) SCA LINDVI 


61-Coefficient Remez Bandpass Filter Design 


Figure 4.3. 
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(DB) 
0.0 


MAGNITUDE 
-70.8 60.0 —50.0 --40.0 -30.0 -20.0 -10.0 


Figure 4.4. 


25664 1.88496 2.51328 
footer Vee tile PA 


21-Coefficient Remez Bandpass Filter Design 
(Band edge spacing 0.01) 
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Figure 4.5 21-Coefficient Remez Bandpass Filter Design 


(Band edge spacing 0.05) 
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Figure 4.6. 21-Coefficient Remez Bandpass Filter Design 


(Band edge spacing 0.1) 
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MAGNITUDE (DB) 


-60.0-50.0 -40.0 -30.0 -—20.0 —10.0 
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0.00000 0.62832 1.25664 1.88496 2.51328 3.14160 
PREQUENG@)] fe iA 


Figure 4.7. 21-Coefficient Remez Bandpass Filter Design 
(Band weighting 1:1: 1) 
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Figure 4.8. 21-Coefficient Remez Bandpass Filter Design 
(Band weighting 10:1: 10) 
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Figure 4.9. 21-Coefficient Remez Bandpass Filter Design 
(Band weighting 10: 10: 10) 
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Finally, Figures 4.7 to 4.9 are the result of varying the weights of the 
original filter specifications as follows for N = 21. 
Figure 4.7- Band 1 weight = 1.0 
Band 2 weight = 1.0 
Band 3. weight = 1.0 


Figure 4.8 - Band 1 weight = 10.0 
Band 2 weight = 1.0 
Band 3 weight = 10.0 


Figure 4.9- Band 1 weight = 10.0 
Band 2. weight = 10.0 
Band 3 weight = 10.0 


Again, as expected, the frequency bands with the heavier weight assigned 
more closely approximated the desired frequency response, Figure 4.8. Here it 
should be noted that the weights are taken into account relative to one another. 
Figures 4.7 and 4.9 are identical because in both cases the relationships between 
the weights is 1:1:1. 

In this example, we have seen that the Remez exchange algorithm provides 
an effective way to design linear phase FIR filters. However, a primary drawback 
to this procedure is that the CPU time required grows quite rapidly (on the order 
of N) with the filter order N. 

For example, fifteen CPU minutes are required for the design of a lowpass 
filter with N = 512 using a CDC 6500 [15]. This is the time required for one 


iteration of the procedure. Typically, more than one iteration is required to find 
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technique to reduce the CPU time required is desirable. Reference 15 presents a 
method whereby the pcoperties of a high-order filter can be extrapolated from a 
lower-order filter. 
B. METHOD FOR THE DESIGN OF HIGH-ORDER LINEAR PHASE 

FIR FILTERS BASED ON A LOW-ORDER PROTOTYPE 

A detailed explanation of this method will not be presented, however, a sum- 
mary of the general idea behind this technique follows. The interested reader is 
directed to Reference 15 for details of its application. It should be noted that the 
term “high-order” refers to filters with orders approaching N = 2048. 

The underlying basis for the high-order filter design technique is the observa- 
tion that in high-order, multi-passband/stopband filters, extremal frequencies (i.e., 
ripple frequencies) in the broad part of these bands are more evenly distributed 


than in the transition regions. 


48 
jaste | 





Figure 4.10. Desired High-Order Filter 


Figure 4.10 represents the desired frequency response of a high-order multi- 
passband/stopband filter; and it is noted that in regions A, C, and E there is a 
uniform distribution of extremal frequencies, while in regions B and D the distri- 


bution is non-uniform. 
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The design of this filter using a lower order prototype involves the following 
steps: 


e To obtain the low-order prototype, the uniform extremal frequency regions, A, 
C and E are “cut-out”. The remaining non-uniform regions B and D are then 


merged to form the lower order filter whose frequency response Hp (e’ *) is 


shown in Figure 4.11. 





Figure 4.11. Low-Order Filter Prototype 


e The Remez exchange algorithm is then used to obtain the filter order, N, and 
the coefficients required that will meet the specifications of Hp (2%). This 


is done by varying the order and the weights, until the order and weights that 
yield the desired passband ripple and stopband attenuation are found. 


e Once the low-order prototype filter 1s obtained, its extremal frequencies are 
plugged back into the appropriate regions of the original desired high-order 
filter, Figure 4.12. 


e The extremal frequencies are used as initial values, and a final run of the 
Remez Exchange Algorithm is then made to obtain the filter order and impulse 
response of the high-order filter design. Again, the filter order and weights are 
varied until an optimum design is obtained. 


In conclusion it should be mentioned that in Reference 21 a program is outlined 
that automatically designs high-order FIR filters using this procedure. Approxima- 


tions of filters with orders up to N = 4096 have been achieved using this program. 
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Figure 4.12. Low-order filter prototype being placed 
back into the appropriate regions of the 


desired high-order filter 


In conclusion it should be mentioned that in Reference 21 a program is outlined 
that automatically designs high-order FIR filters using this procedure. Approxima- 
tions of filters with orders up to N = 4096 have been achieved using this program. 

This section has illustrated the usefulness and relative ease of application of 
Computer-Aided-Design (CAD) techniques, such as the Remez exchange algorithm. 

Since this algorithm is only applicable to linear phase FIR filters, it is worth- 
while at this point to diverge to a discussion of a popular algorithm for the design 


of recursive IIR filters, the Minimum p-Error Design Method. 
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C. THE MINIMUM p-ERROR DESIGN METHOD 

Before proceeding with a discussion of this technique, it should be mentioned 
that, in the interest of brevity, the pertinent equations and relationships used 
are presented without elaboration. This is because the mathmatical basis for the 
minimum p-error criterion is quite extensive and beyond the scope of this thesis. 
The interested reader may find details of its derivation in Reference 19. 

The minimum p-error design technique consists of a generalization of the min- 
imum mean-squared error design method, wherein the optimum filter coefficients 
are determined by minimizing the following error function [18]. 

k 
E,(A) = ) | w (9x) [H(A, x) — Ho(4x))” (4.9) 
k=1 
where 
w (6,) = weighting factor 
H(A,6,) = approximating function (actual response) 
Hp (9) = desired frequency response 
6, = digital frequencies over the range of interest, k 
A = vector containing the k independent parameters 


(i.e., the filter coefficients) 


Note: If the value for p in Equation (4.9) is unity, this relationship describes the 
mean-square error. 

The filter approximation problem can thus be stated as follows: Given an 
amount of error that can be tolerated, E2,, find the set of k parameters (filter 
coefficients) such that the error between the approximate frequency response and 


the desired frequency response is within the stated tolerance. 
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Again, looking at Equation (4.9) it can be seen that for large values of p, 
the approximate frequency response, H(A, 6,), has to be very close to the desired 
frequency response, Hp(6;), to remain within the pre-selected error tolerance value 
of E2,. Furthermore, a sufficiently large value of p results in an optimal solution 
that is very close to the optimal Chebyshev (or minimax) solution [19]. The method 
of solution of this approximation problem, as will be shown, is dependent on the 
error tolerance value selected, E2,, the form of the transfer function, H(z), and 
the parameter vector, A. First a discussion of the transfer function is presented 
because its form has a direct impact on several important features in digital filter 
design [18]. 

As is well known, three possible forms of the filter transfer function are direct, 
cascade and parallel realizations. With respect to errors directly related to the 
physical structure of the filter, i.e., quantization effects and finite coefficient size, 
the use of the direct form of the transfer function in filter realization is undesirable. 
This leaves the parallel and cascade forms, with the cascade form being selected 
because the zeros of the transfer function remain unaltered in the process of cas- 
cading, resulting in well defined stopbands. Conversion of a filter transfer function 
to the parallel form, on the other hand, involves partial fraction expansion, which 
results in the zeros being less well defined. Theoretically this should not be the 
case, but, when the parallel form is implemented digitally, the zeros are slightly 
different than those of the original direct form transfer function due to the finite 
precision of the computer. 

Using the cascade form of the filter transfer function has other advantages 
including: 


e stability tests of the filter can be accomplished without lengthy calculations, 
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e the frequency response is of a simple functional form that readily lends itself 
to insights as to how the poles and zeros (hence, the filter coefficients) impact 
the error function. 


Thus, the first step in the use of this design procedure involves decomposing 
the proposed approximating transfer function, H(z), into cascaded second-order 


sections. 


H(z) = bo |] Hil2) = bo T] ee 4.10 
(z)= off (2) =*o [T S55 ah. (4.10) 
=1 t=] : : 


where - _ 
ko = a positive normalizing constant, and 


N = filter order /2. 
The parameter vector, A, therefore, consists of the multipliers used in the 


cascade structure. 
A= [a11, @12, by o12;- - +5413, 423, Din ba;, CO Oy ko] (4.11) 


The following expressions for the frequency response and group delay of the 


filter incorporate the cascade structure and parameter vector. 


e Frequency Response 


N ; 
1+ a2;)cos9 + a3; + 7(1 — a2;)sin?@ 
PUNE 8) = Palaces ectenatra7s)smF 64 15 
ey, te") oT] (1 + bo;) cos + b1; + 7(1 — b2;) sind ee) 
e Group Delay 
N Pa: 
d 2cos 8+ b1; + j2sin 6 
A 6 — j@ = po IE —E 
Ts Se ry sla Tals co Segerer eee 
sige 2cos 8 + ai + j2sing | ) (4.13) 
(1 + a2;)cos? + ay; + 7(1 — agi) sin 8 


Equations (4.12) and (4.13) describe the frequency response and group delay 


of the approximating function, H(A,6,). To determine if the coefficients of the 
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A vector yield an optimum filter design based on the desired frequency response 
Hp(6,) and the error that will be tolerated E2,(A) involves finding the partial 
derivatives of the frequency response and group delay portions of the error function, 
in terms of the filter coefficients, and minimizing them. 

In other words, the filter coefficients comprising the A vector are to be selected 


so as to minimize the following partial derivatives: 





ae = yp (64) = oe — (a (A, 9x) lp (8,))°?-* (4.14) 
Jean Yop () = =G (A,6.)-— Hp (ene (4.15) 


and similarly for 0F2,a(A)/06,; and OE2)7(A)/0bj;, ete. 

Upon examining Equations (4.12) and (4.13), it is that apparent these par- 
tial derivatives are not easily attainable because complex expressions are involved. 
To remedy this the polar coordinates of the poles and zeros are used for param- 
eters, rather than the original rectangular form. Thus, the parameter vector and 
expressions for the amplitude and group delay are modified as follows, 


e A Parameter 
A= [ro ’ 901 ’ Tals Upls 220 9T0i5 90: lpiy One bes 9 ko| (4.16) 


e Frequency Response 


aA, 8) = ke py sense eRe inel iat eoth ata pene) 


{1 — 2rp; cos (6 — 65;) + rpi2 [He {1 — 2rp; cos (8 + Opi) + Tpi2 
(4.17) 
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a 


e Group Delay 


N 
FA. 6) = Ss | 1 —r,; cos (@ — 6,;) 1 —r,; cos (8 + §,,;) 
i=! 
ms 1 — r,; cos (8 — 4,;) 1 —r,; cos (@ + 4,;) 
1 — 2r,; cos(6 — O4;) +ro;2 «1 — 2r,; cos (8 + 95) + To x2 


resulting in the following partial derivatives of the error function: 


OE,,7(A) < Or ie 
Na a oan p—1 
ar. = 2d Pwr (6, ) a (r(A, 6.) H p(,)) 


and for 0E2, a(A)/06,; and OE, 7(A)/069;, etc. where 


Ca E Tosa COS (6 a Oo) 
Oro: me 2r oi cos (8 — O5:) + Toi? 
1 — 2r,; cos (6 + 6o;) + ro:2 
Oo _ ,|___raisin(9— Ooi) 
08,; = 1 — 2r,; cos (6 = OG;) + 1 o2 
Toi sin (6 + 4; ) 
1 — 2r,; cos (6 + 695) + rox2 
OF ~ (1 et Tpi2) cos (@ — Ga.) — 2rpi 
Or pi (1 — 2r,; cos (6 — 6,;) + rpiz) 
(1 — yi) cos (6 + 65;) — 2rp; 
(1 — 2r,; cos (6 + 65:) + rp) 
Or = Ti (1 — T yi) sin (6 = Ont) 
oi (1 — 2rpicos (9 — pi) + Tpit)” 
T pi (1 — Tyi2) sin (0 + Oa) 
(1 — 2rp; cos(6 + 6,;) + pz) | 


os 








1 — 2r,; cos (0 — 85;) +'Tp:2 34 1 — 2r,; cos (0 + 8,;) + rpi2 


(4.18) 


(4.19) 


(4.20) 


The partial derivatives 0a/Orp;,0a/00,;,07T/Oro;, and Or/06,;, are the same as 


the above but with changed signs. 


Thus, the approximation problem amounts to the minimization of nonlinear 


functions (Equations (4.19) and (4.20)), of n variables (the poles and zeros). This 
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problem is readily solved through the use of the Fletcher-Powell algorithm, the 
details of which can be found in Reference 20. 

A program that performs the synthesis of recursive digital filters using the 
minimum p-error criterion and the method of Fletcher and Powell for function 
minimization is contained in Reference 16. Included in this reference are an exten- 
sive program description, input requirements, dimension restrictions, and examples 
to illustrate how the program is used for various types of filter design. Also included 
is a copy of the actual program. 

The reader is reminded, however, that to make use of the program, the transfer 
function must be in the cascade form mentioned previously. The program output 
consists of an argument vector, X, and a gradient vector, G, corresponding to the 
minimization of the functions Fl, F2, and F3 which are used by the program for 
the magnitude approximation, the group delay approximation, and the combined 
magnitude and group delay approximation, respectively. Based on X, the poles, 
zeros and coefficients of the cascade realization of the filter are computed. If desired 


by the user, the frequency response is also given [16]. 


D. SUMMARY 

In summary then it is appropriate to cite a recent paper by Little and Gowdy 
[17], that evaluates both the optimum FIR design method (that employs the Remez 
exchange algorithm), and the minimum p-error method used in IIR design. Their 
evaluation specifically investigates convergence problems encountered when using 
these iterative techniques, and considers the design of lowpass, highpass, bandpass 


and bandstop filters. 
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The following is a brief summary of the more important points mentioned in 
this article, that should be considered when using these iterative techniques. 


1. Minimum p-Error Design Method 


a. Advantages: 


e Extremely flexible - can be used to design filters with arbitrary magni- 
tude and/or phase characteristics, as opposed to being limited to lowpass, 
highpass, bandpass, or bandstop designs. 


b. Disadvantages: 
e Nonconvergence problems due to: 
o a poor guess for the initial parameter vector X 


o finite wordlength 


e Uses large amounts of CPU (IBM 3081) time (up to 2 minutes for higher 
order filters, greater than N = 12) 


e Requires a considerable amount of input to specify one filter 
2. Optimum FIR Filter Design Program 
a. Advantages: 


e The program is easy to use 


e Consumes relatively little CPU (IBM 3081) time, provided the filter is 
not of extremely high order. 


b. Disadvantages: 


e Restricted to the design of the more common filter types: multi-band, 
bandpass, Hilbert transform, and differentiators 


e Highpass and bandpass filter designs were poor when even filter lengths 
were used, but good for odd filter lengths 


e The user must supply an FFT program to obtain frequency response data 
to verify the filter design. Convergence of the program does not guarantee 
an acceptable design. 
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Keeping these considerations in mind, both programs can be used effectively 
to design a wide variety of FIR and IIR filters. It is suggested that if a designer 
wishes to use either of these programs this paper is a valuable reference to assist 


in identifying problems that may be encountered. 
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V. CONCLUSION 


The areas of recursive (IIR) and nonrecursive (FIR) filter design have been 
investigated, while the more predominant design methods have been extensively 
discussed and exemplified. Additionally, the prevailing computer-aided-design al- 
gorithms (Remez exchange and Fletcher-Powell) were also presented. 

In conclusion, points in these three areas that merit special emphasis are 


summarized below: 


A. RECURSIVE FILTER DESIGN 


e The desired frequency response may be obtained using a lower order filter 
than if a nonrecursive realization were used, however, filter stability must 
be considered, and a linear phase characteristic is not guaranteed. 


e The traditional design methods involving Butterworth, Chebyshev or el- 
liptic analog prototypes, and the bilinear transformation are algebraically 
intensive. The direct design method presented eliminates the need for 
determining an analog prototype and for prewarping, thereby reducing 
the number of calculations required. 


B. NONRECURSIVE FILTER DESIGN 


e Although a higher order filter is required than in recursive realizations to 
obtain the same frequency response; nonrecursive filters are always stable 
due to the all zero nature of their transfer functions, and can exhibit a 
linear phase characteristic. 


e The filter coefficients, h(n), may be determined analytically by expanding 
the desired frequency response in a Fourier series, because the Fourier 
coefficients correspond to the filter coefficients. 


e Windows may be used to eliminate the overshoot in the frequency re- 
sponse (Gibbs’ phenomenon) caused by truncating the number of terms 
in the Fourier series to a finite value. However, windowing decreases the 
sharpness of the filter’s cutoff region. 
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e Frequency sampling can be used when an analytical expression for the 
desired frequency response cannot be found. The IDFT is then used to 
determine the filter’s impulse response from the these sample values. 


C. COMPUTER-AIDED DESIGN 


e For FIR filters, programs using the Remez Exchange algorithm are the 
most popular, while for IIR filters, the Fletcher-Powell algorithm is in 
predominate use. 


e CAD is especially advantageous in the design of extremely high-order 
filters, or filters with arbitrary frequency response characteristics. 


As stated in the introduction, the design methods presented here are by no 
means the only ones available to the filter design engineer; however, due .o the 
quantity of information available, one is easily overwhelmed. A need exists for a 
concise source of information on the popular design methods available, and how 


they are used. This has been the intent of this thesis. 
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APPENDIX 





MAIN PROGRAM 


MAIN FROGEAMs FIFE LINEAR PHASE FALTER CESIGN FRCGSAS 
ADTHORS: JAMES H. SCCLELLAN 
DEPARISENY OF ELECTEICAL ENGIUNEEFING END CONPULEF SCIENCE 
BASSACUUSETIS TUSIZTUSE Ce TECENCICoY 
CAMBRIDGE, MASS. 02135 
THOZAS W. PARKS 
DEPAKIMENT OF ELECTRICAL ENGINEERING 
BICrE JNIVEESTIITY 
HOUSTON, TEXAS 77001 
LAWRENCE R. KABINER 
BELL LABCRATCRIES 
MORBAY HILL, NEW JEBSEY 075784 
INPOtL: 
NES Pde ier LENG Le 
Dliro== (YPE Crerceier 
1 = MOLTIPLE PASSEAND/STOPBABC FILTER 
2 = DIFFERENIIaTCF 
3 = HILBERT ITRANSFCReH FILTER 
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Marin=— GalD DENSITY, alll EE SET IC l6 UNLESS 
SPECIFIED CIHERWISE BY A POSISZVE CONSTANT. : 
EDGE (2*NBANDS) -- BANDEDGE AZRAY, LOWER AND OFPER EDGES FCF EACH BANC 
SITh 4 MAXIHCY GF 10 EANDS. 
EX (NBANDS)=— CESZ&ZED FUNCTICN AsEAY (Cz PESIRED SLCre IF A 
DIFFERENTIAICR) FCc EACE LEADL. 
WIX (NDANDS)-—- @ELGHT FONCTION ARRAY IN cACH BAND. FCS A 
DISS EFENTIATGR, I65E wEiGEL FUNCITICnw 15 INVErSELY 
PROFOaTIOWAL iC F. 
SAMPLE INPUT DATA SESZUP: 
aap leo, ul 
Pavel, 02 0035 
(6645 5,055 
Noe meee 
TOSU SO 10s 0 
Topo pha oPecisiEsS & LeNGik 32 BANDPASS FILTEL @128 
STCPHANDS 0 TC 0.1 AND C.425 TC 0.5, AND PASSBAND FRCS 
Jec 40 U.35 sITH WEIGHTING OCF 10 Js IHE STOPBAsDS AAD 1 
IN TBE PASSBAND. THE GPID DENSITY DEFLULIS IC ilo. 
Toes 2> DHE FILIER IN FIGURE 0. 
THE FOLLOWING INPUT DATA SEFECIFIES A LENGIH 32 FOLLEAND 
DIFPFERENTIAICR &8ITH SLOFE 1 AND WEIGHTING CF 1/F. 
THE GRID DENSITY WILL EE SEI TC 20. 
32 g2,1,20 
Oye 
1.0 
1.0 
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eAD,DEV,1,Y,GEID,DES,8T, ALPHA, EXT, NFCNS, NGEID 
PS/NITES,LCUT 
7 AD (66) , ALEHA(66) ,X (66) , ¥ (66) 


D(1085) ,wT (1985 
0), 81X (10) ,DEVI 


Alkp >! 


kr (10) 
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B, 182, 1BN, 1HD/ 


On_ON0NNtyouRH 
mAs fib ICI tH th 
—~* MmrUnMNmnHMc td 
atWtIkPMEtK. od 
~—“cncot.o™., 


167 


2) 


AANANAANNA 


Anan 


NANA 


QANNDA 


QANANN 


ICUT=11MACH (2) 
P1=4,0*DATAN (1.000) 

Fiz=z-OD00*P1 
TBE PERCGRAM IS SET OP FOE A MAXIMUM LEKGTIE CF 128, BUT 
THIS UPEEH LIZI2 CAN BE CHANGED £Y BELISENSICNING THE 
ARAAYS iEXT, AD, BLPHA, X, Y, H_TO BE NEMAZ/2 + 2. 
THE ABAAYS DES, GRIL, AND WI BUST DISENSICNED 
lo (NFMAXs2 + 2). 

NPMAX=128 
100 CONTINUE 

JTYEE=0 


PROGRAE INPUT SECTICH 


READ (4,110) NFILT, JTYSE, NEANDS, LGRID 
IF (NEIET. 20.9) SICE 

110 POBAAT (4 (15,3 } 
IF (MPILI.LESNPSAX.AND.NFILI.GE.3) GO TC 115 
CRUEL OERECE 


SICP 
115 IF(NBANDS.LE.O) NEANDS=1 


GRID DENSITY Is ASSOZED TG BE 10 UhLESS Srecirier 
CTHEHWISE 


IF (LGRID.LE.Q) LGRID=10 
JB=24hBANDS 
REAL eo CF a SL NS AHO 
BEAD (4,120 FX {J) ,J=1,N3ANOS) 
READ(4, 120) [Wt (J) J=1,4BANDS) 
1Z0 FOGMAT (4F15.9) 
IF (JIYPE.GT.O.ANL.JTYSE.LE.3) GO TO 125 
CAal ERKCE 
SUSTGE 
125 NzEG=1 
It (JTYFL.~EG.1) NEG=0 
NCLD=ENFILT/s2 
nwODD=NFILTI-2*NCOD 
NFCKS=hFILT/2 
iF (NCDD.EG.1.AND.NEG.EC.0) NFCNS=NFCNS+1 
SET OP THE LENSE Gall. ZHE NOMESR OP POINTS SIN TERE GREED 
Is {FILIEX LeENead + 1) "S620 Gensa7272 
GRID (1) =zDGE (1) 
DELFHLGELDENFCNS 
DELF=0. 3/DEuE _ 
LE (NEG22y.0) SOIC, 135 
‘ TF (EDGE (1). T-DELF) sRID(1)=CELE 
135 CONTINUE 
J=1 
1=1 
LBANL=1 
14Q FUR=EDGE(1+1) 
145 TEMP=GE1D (J) 
CALCULATE THE DESIREL MaGNITUDE BESPCNSE AND THE WEIGHT 
FUNCTICN CN THE GhIL 
BT (OOM RAE (Le we Ea eee LOA ae ERE 
ei eiteee LEME, FX,artX,LBAAL,JTYPE 
=J+ 
GID (J) “TEMP *DELE 
iF (GRID (J)-GI.FOP) GO IC 150 
GC 1C 145 
150 GKID(J-1) =FUP 
Bee eae kee Ue ry eae TEAnT UTEDE, 
a a esl te kee ite haps) b er 
LEEND=LBAND*1 
=ltegz 
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